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ABSTRACT 
The four estimation algorithms in the University of 
Alberta Identification Package were tested using 
computer generated and experimental data. A double pipe 
countercurrent liquid-liquid heat exchanger and a binary 
distillation column with eight bubble cap trays were 
the two experimental units used in this study. An IBM 
1800 digital computer, executing Direct Digital Control 
(DDC), was used to control the units and to accumulate 
the data. 
A considerable portion of this thesis is concerned 
with the theory of the four identification routines: 
1) Generalized Least Squares (STIME) 
2) Maximum Likelihood (ASTRM) 
3) Instrumental Variables (IVOFF) 
4) Unbiased Estimator (UBIAS) 
and the software of the package. The University of 
Alberta's Amdahl 470V/6 digital computer was used to 
implement the scheme. 
The résubts ote’ eth e computer generated data 
identification runs demonstrate the effect of: 
1) lengthening the input/output data sequences 
on the accuracy of the parameter estimates 
2) increasing the model order on the accuracy of 
the parameter estimates 


3) the type of identification scheme used 
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4) the effect of noisy data on the accuracy of 
the parameter estimates 

These results were compared with findings from other 
studies and general agreement was found. The 
generalized least squares identification routine 
produced results that were unacceptable for good process 
modelling. The preferred identification routines were 
the instrumental variables and the unbiased estimator 
methods. 

The experimental results obtained from the heat 
exchanger were compared to previous investigations and 
the models estimated for the distillation column were 
used by a fellow student. Some problem occured in 
estimating the gain of a process with a large time 
Constant. A solution was proposed and proved with 
testing, that this problem could be ‘reduced if the ‘clock 
interval of the input disturbance signal was increased. 
The model proposed for the heat exchanger was in good 
agreement with models from other investigations but some 
problems were encountered in estimating the time delays 


of the process. 


cw het 


Lae “ie Loot, 
o e s a? tae 

7 i i” yA 

» 

hte : 


Yo. yoexuges ans tO” eteb eter Fe soto 4 o. 
cise tt Sen 7e79m8 7189 oe = a 


~— 

























| ad 
sense vet? arbiint) .iie-w bs Tegmas 910" éstyaer- e 


eat .bnuct’ ee” jnenseaee feraneg sae es] 


. 


“ : ¥ 
~ S's 


\ C : ; 
sujenn Anelredéeitnet., (ee7aeps juseol hesilersneg 


ekegerd, poaog It Sidedere~ sPy sem 26d? OF i 4so% bea 


wAluvengt *hiawebi pews ts 9 @ut ~gndst 


giew 2an?qyor 


ssgamprise veecetdry, gfe tne gaicoi my fasacwungens 
+ 80d 


ejjivean Letaemiveqss sat 


oo. iw bs e79 
Swe enol iegbtac ne solvate) ot beipadco ete¥ vege 

syvew oaYos p»AeLebd Cake Riy Arges cade aie eisboms 
wv! wont “aeideng ea jnebugs woftet? « Xe 

niv? ez a “det eescorve Bete cqica sit porseahde 
dere bevona #08” Be2cqgotg ~ Sey fetdpy as A _.dneda 
Jouyig 684. tk Seb hes ec »sfuod wattend afd? teas ,aahde 

‘ 

besesvual caw Zenyle aonadyeur te juget ¢e2% Yo ley 


row “eAdedoxe Jest ens uot tareqntq fabom 


stddevet rote wot S2ekta Wily os 


enae ud eeoetssEancs 
- : 4 5 i 
eyaish sit ohe Recser. ‘se at Us sJ7og0005 stew 2° 


: . .eee0o1g 2 ad 





TABLE OF CONTENTS 


CHAPTER PAGE 
a gs Jigaleaveae relieves. Geile)” eS ees. ot Sao eo ae 1 
ois Theory of Identification Methods 4 
MIC ROGWUGC DON =... 6.54. «6 be He Fe Me 4 

Statement of the Problem 5 

PROUT OW UTUS HO Ce Bere te Be He te H's g 

Generalized Least Squares... 10 

MaxmwmomneuL kel inood “e. «. «. b. ++ « 21 

Divs ume calwVarLaDles  . 20. 30 
Unppbascdakstimatorfete.s.for. . a3 

a Mee Organization and Structure of the 


University of Alberta Identification 


es el) ao oe Eg ae Bee oe 0 38 
Life \Sefede thvers 5 elon; Ja ix oe i are ce ee 33 
BeOrnameMOdul@s@AGuire . «. +. +. « 42 

ODEs Gs mee. VOLUN w+. 0. e. 42 
Dink Hoesen BOS. 8. ae He ee HS ee 48 
PO Ls PARE cele cams es B+ #5 ioe oo e's 56 
DATGE VRP AOR. COP UMD 2 ke ee os 60 
DAGEN sous) 2) eee eee ee eee 64 
DMO His @camic is e's A om sys 400 6° Bis 74 
CHCOnMME EE tie tows Do ngs Sr 84 
SD PIMERSSOUMAVIALDOOR 2s. +. 86 
DS Meek 54 ee.) as oe, a8 be ai le 97 
IVOFF oe Oe ee Pe eee 106 


UD DASE tee sete em oth a ee ve! tO 
DIAGN ee eee 4c ee eS wes te te AEE 


MD US Fae eee sere fo 7S led 
NO. LOG wtn UR phils Been e's be Ue ke 9 
Bu NA lisse ee”. ce a 4 


iv 



















. enrtoubesie® . 


if . - . - « * . . . . 


Ps 


( 
B + \\ebantemM doltpettianess te yrosst + 


bo ae oe nagzoebor2 Tl 


gaicoit eat Ie freare to 72 im 
c gd ae ee Ae spol (40a : 


6 . 2 ees bpe, fepes fexiisteacd. 
hoor [adi goede 


= ois « . * . ’ 
O€ ; . eels. se in tvewmpers eal 
-F 5 Sa nojeurtsed besedtdot. 


aft 1d srugev 7c tre wotersataegatO . 
ao tgeeiasianact atrediw 6 vyoiereavinet | 
4. shadond 

aci@gubotsal 

=i . , «9° fete bos Gs55o7%% 
7. fis S ‘Sl ee eae A 
by! fe Sp < he ee ee 
jc “ae ee pe): Lone ee .- ae 
he : ¢. ocen a Oe eee 
i Th aApa ge ayamd 
| ot, be A. eee eee 
gi’ (2. BE ee See 
“hiZe 


CHAPTER 


is Bes 


Ve 


Wook fe. 
CU See ee 


Parameter Estimation of a Simulated 
Process using UNALIP 

Introduction 

Test Cases 

Identification 

Summarized Results 


Discussion 


Parameter Estimation of Models for 
Perot ance Units 
Introduction 
Distillation Column 
Heat Exchanger 
Experimental Procedure 
Distillation Column 
Heat Exchanger 
Results 
Distillation Column 
Heat Exchanger 


Discussion 

eeratae co send Recommendations 
Bibliography 

Appendix A 

Appendix B 


MOOG Gee ees > 4 ee eG 


PAGE 


rae) 
145 


150 
Nyse) 
ioe] 
155 
is3 
176 


178 
178 
GS 
181 
183 
183 
183 
185 
1O5 
194 
ee 


mic 


206 


216 


eee 


fdmiz se to noktens tea 
41JAWY yoter seec07% 





-_- «a ee yaetandat 
—o7 appe oeet 
noldaniyaageo 
‘(yeu oo gd tenes 


ADE FL. so fangoesd 


, » 'pobhspty rage 
age Fs ) no bee fi tiene 


amuieo ooldgesg twa tI 


. tentetoe’ ert 

















yoseaeaat 


oM te ‘petitenige® Seoteussst 


fia fostd we1is 


a « aoe « wepeamoeeta 





aclisdo spore bos ecobewionss 


JeWoeK=Fefe hiya 10); : is a ar ~ eh2 


vi 





LIST OF TABLES 


High Noise/Signal Ratios for each 
Noise Sequence and Process 


Low Noise/Signal Ratios for each 
Noise Sequence and Process 


Comparison of the Standard Deviations 
of the Error Measures for Generalized 


Least Squares (STIME) 


Comparison of the Standard Deviations 


of the Error Measures for Maximum 
Likelihood (ASTRM) 


Comparison of the Standard Deviations 
of the Error Measures for Instrumental 


Variables (IVOFF) 


Comparison of the Standard Deviations 


of the Error Measures for Unbiased 
EStumatore CUBLAS) 


Comparison of the Average Standard 
Deviations of the Parameters 
Estimated by Generalized Least 
Squares (STIME) 


Comparison of the Average Standard 
Deviations of the Parameters 
Estimated by Maximum Likelihood 
(ASTRM) 


Comparison of the Average Standard 
Deviations of the Parameters 
Estimated by Instrumental Variables 
(IVOFF ) 


Comparison of the Average Standard 
Deviations of the Parameters 
Estimated by Unbiased Estimator 
(UBIAS ) 


Comparison of the Computational Time 
(seconds) Required for the 
Estimation Routines 


PAGE 


160 


160 


163 


164 


165 


166 


167 


168 


169 


170 


Ailes 
























aids tc? sot aa tpeute actos A 
wai ca923~%9 bre sogeupse S% 


\ ¢2g8 301 tots Ba fangi@\es rok, Woe 
rer soanodtd Boe aoneupes otis 


re : + . 

sHobiel ves Pashia de ond Yo Ritts plage * 

‘ bee S + yatey “ot 25 erent WY -O7F* gy sag @ 
[anit Regs) A jeasd 

pnoige bys trebesrt? (etd to meekrigaod 
ure ted tok fon of yours od Io. 
‘RTSA) voo6ood rieail 


phe dz sat to noetteqaed 


2, an 
7 
~~ hy, ae ¥ ve ran TOS ertis to 7 
aot . (350Vi) aefdsite? 
avghdts tryed aye it Yo netivegmoe? 
Si “> & how ad >i’ J 547 2 b o > 4 
: mi : . 
Ss 77. 28.9 sM “oyvyS FAs FO 
(2al Sse? sosentses 
’ 
ovebrgJge spaitevk 2hI80 poet yaaees 


4 


>t epye® shy Vo sateteived . 
of &.] : $2e2- 44. 2 denises 
ror. _CUSMIT3* <9 yeupe 


4 
¢ 
a 

5 
» 


piphme 22 eRe rava =4/ o-ont tageee? 
tvademe se 3 tos Fo anodes gad 


Nona treads: ng keh. 48 nee 

. * CMaTSA) 5 

i“ebiiete® Spe 1gva vt? To set sipqe@od 
Amster ‘he. Yo. Stt902 vee: 

aStigmuA teri ¥d bed rapitet 
+. GORE) 


yeone? 


tvs@ 


“n 

a 

ry 

wh 

7. 
: 

mh 


 dvabts72 Sap ye vd a2 Io Moe VEG 
oyedatete? acxv % plekeled 
Lad acl ident =. S636 







TABLE 


Comparison of Parameters Estimated 
by Three Identification Routines for 
the Steam and Bottom Composition 
Transfer Function 


Comparison of Parameters Estimated 
by Three Identification Routines for 
the Steam and Top Composition 
Transfer Function 


Comparison of Parameters Estimated 
by Three Identification Routines for 
the Reflux and Bottom Composition 
Transfer Function 


Comparison of Parameters Estimated 
DY ihree Identification Routines for 
the Reflux and Top Composition 
Transfer Function 


Identified Gains of Simulated System 
Using Each Estimation Technique 


Effect of Clock Interval on the 
Estimated Gain of a Simulated System 


Vi. 


PAGE 


186 


188 


189 


191 


ied 


ogT 


- 





is Vere 
‘ae DP 
ae | 
‘bevemtdes tindvue el to noses 
ao) asrhjyew sg: piawael sem 
Ao bsLE Sg stod © bee mee 


‘Dedetehded roa Yo scott a. 

192 gealtvont ne¥eel i gsnebt sath? ga 
coistteagaod. gel voa meege © 
AOltonws 74 e084 


j 
badamtstd ec s}etiersi§ jo apes yer 
"Ot @yat-voel agteesliizgead( sac8 “ed 
ratibergae? yyssod Fad, sulted bets 

| norjchet yetadedt 


betemist=3 ec tome i89. to cgeetaeqeed 
10% gonigvc® wobjestlbsass! sssat 9a 
nolglecgmt? aok fae sy ited png 
eiigened aetayeet 


Civey? Setgfoaes te 2afew Hebei Snedt 
suprofost neoltanlied does antel 


siz oo ([sveesal aAocths th dost 
metzyc osisdimit @ to ptap betenkded 


ee 



















is 5 a 
—) 4 


es 
— i. 


Sc 


ss 








= 





FIGURE 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


LIST OF FIGURES 


Diagram 
Diagram 
Diagram 
Diagram 
Diagram 


Diagram 


for 


for 


for 


for 


for 


for 


twee 


STEER 


STEER 


STEER 


HELP 


DATGET 


Box Diagram for DAGEN 


Flow Diagram for DAGEN 


Flow Diagram for DAGEN 


Box Diagram for DMOVE 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Flow 


Diagram 
Diagram 
Diagram 
Diagram 
Diagram 
Diagram 


Diagram 


Diagram 


Diagram 
Diagram 
Diagram 
Diagram 
Diagram 
Diagram 


Diagram 


for 


for 


for 


for 


for 


for 


for 


for 


ror 


for 


for 


for 


for 


for 


for 


DMOVE 


CRCOR 


STIME 


STIME 


STIME 


ASTRM 


ASTRM 


ASTRM 


IVOFF 


IVORE 


UBIAS 


UBIAS 


UBIAS 


DIAGN 


MDISP 


x 


PAGE 
3 
49 
50 
5 1 
5"t 
61 
66 
67 
68 
(ps) 
76 
82 
88 
89 
90 
ag 
100 
101 
107 


108 


119 


126 


c 
Fé 
t 
{ 4 
i« 
4% 
ra 

~ 

J 
i4% 

1 
Fes 


oy Poa 
Wig % 
— © 


hele. Zi 


TIH6 Act 
HaSTS et 
s2HCe 02 


io 


AASTS tO? 


TuaH at: 


TRDTACG tcl enon [a wot’ : 
i 
SIA s01, mayeetd tee + 
340 
taOA io% Siserd woht 
a 
AA a p 
aa i ror pe eae Ae & cee ae 
~ 
iC 1 
SYVUMY FOL NEN RE ES, Bos : 2 
i : A 
51 wetwerd wei 
¥ A “az ws" LBS. iT wold 
ci 4 pet hake woL? ' 


VATEA 02) 
MATA) to3 
HERE 9 
4 TIVE 102 


> 
TIOVI, +02 mabe se wort a 


GA PHS x0 terete a ial 


Sarau ” 
alae 


i 
—matge ra wold 




















Ai ‘ 


7 
=e a g) 
7. y 





; 
; * t 
nanasit ist a 7 
ss a: 
e 


’ ay : ~ 
Peeey ea voit " 
asigeih base 


ph HELD wal 


nA FAs} 4 wolt ; 
Bias ud wor a 
cs Ygsid wold ar 
a Ps 

we anid welt 


geige id wort 5 


875 * ri id 
; ine 


i - 


Siiead: ef 


i letinlis ™ 





FIGURE 
3-26 
3-27 
3-28 


Flow Diagram for NOISE 
Flow Diagram for FINAL 
Flow Diagram for FDISP 
Flow Diagram for FDIS2 


General Block Diagram of 
Simulated Processes 


Pseudo Random Binary Sequence Used 
as Input for Simulation Runs 


Simple Schematic of the Distillation 
Column 


Simple Schematic of the Heat 
Exchanger 


Tube Flow Rate Response for Test A 
Tube Flow Rate Response for Test B 


Shell Outlet Temperature Response 
for Test A 


Shell Outlet Temperature Response 
for? Tests 


Reflux Flow Rate Response for PRBS 
Disturbance of Length 63 


Reflux Flow Rate Response for PRBS 
Disturbance of Length Double 31 


Steam Flow Rate Response for PRBS 
Disturbance of Length 63 


Steam Flow Rate Response for PRBS 
Disturbance of Length Double 31 


Bottom Composition Response for 
PRBS Disturbance of Length 63 
using Steam 


Bottom Composition Response for 
PRBS Disturbance of Length 63 
using Reflux 


PAGE 
134 
183 
140 


146 


1533 


156 


180 


182 
218 


219 


220 


Zen 


224 


cao 


226 


ohh 


228 


eed 


7 ett “ > 9 eee el 
aan | : j - 


sa a ph 
2085 7 ar, he 
beg S210" «ct as 0send: ad oe 



















eff JAWPI 263 es eb wok 


oD 219% Ad} geegorde oe 


( 
; SZid3 a9t eangsid 208% 7 
me 7ipatd vido lt fexvaded 
hy ms ? oy b 4 _. , ' 
a i G41 Zaeuo > [eee Tome ; 
cf sopeuvpeesvaseclis sobtneh obweat _. 
r ecm qofeelval’ Yeh IJwent Ge 
t4ilttalC eds o vitePenie sigarke 
Amy [69 
jeat. ony. Tor $itewene> lease b 
7 -* qegeeroze 
A %@5f ago sznoqash ote walt bdgT 
7 oy . a — > age2H 20 (9 agaT 
pBase i 3 "squat gsijwO Lisége 
A jfeT rot 
A 
2% sJutens(i dot tedzo@ ttsag 
“3 ‘feet’ 1604 
-S43 ol senequel wagh ded) £ekish 
HS 58 ooane! To soocddbaetg | 
oHad 19? ae i stph wold selien 
eS lé slcacd tage ic sugasm debe . 
~ — ® Vv * 
92 : - : te 
CSR 15% ea 2i54 word ar 
a5 (@agued 9 eo: asngveta : 


wie is 
CRAS ot oonadt eof 9 tf * hes 
rssh TE 6idoo’, doga50 bid hatin te 


rs evo 
oe , 





\ note ana ome acaba Fe) 
77 sfpaeis te =o oa 






¥ aot 
2 ee 
es 


4 


es 


| ‘ “ to% Senogeni, re Lada 
¢3 tdemey teu . 










FIGURE 


C-16 


Bottom Composition Response for 
PRBS Disturbance of Length Double 
31 using Steam 


Bottom Composition Response for 
PRBS Disturbance of Length Double 
31 using Reflux 


Top Composition Response for PRBS 
Disturbance of Length 63 using 
Steam 


Top Composition Response for PRBS 
Disturbance of Length 63 using 
Reflux 


Top Composition Response for PRBS 
Disturbance of Length Double 31 
using Steam 


Top Composition Response for PRBS 
Disturbance of Length Double 31 
using Reflux 


Reflux Flow Rate Response for 
PRBS Disturbance of Length 127 


Steam Flow Rate Response for PRBS 
Disturbance of Length 127 


Bottom Composition Response for 
PRBS Disturbance of Length 127 
using Steam 


Bottom Composition Response for 
PRBS Disturbance of Length 127 
using Reflux 


Top Composition Response for PRBS 
Disturbance of Length 127 using 
Steam 


Top Composition Response for PRBS 


Disturbance of Length 127 using 
Reflux 
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CHAPTER ONE 


INTRODUCTION 


UVeE=aLOe Une —"rapiue £rowth’ sor “process ~ control in 
mocerne industry ;**tnere”~ has”*-been” a greater effort to 
increase the ease and the accuracy of model 
Poeentwricatrion™ yincema process “model is Wwsualbly 
required before most of the modern control schemes can 
Demapprred, process identification “is very important to 
PONCLOL, istravewies and Small errors in the model can 
lead to large control problems. Therefore good model 
estimates from noisy data is a must in modern industry. 

This thesis deals with the theory, the structure 
Simerivemmarpr Cavrone Of ea "group or “identification 
programs known as the UNiversity of ALberta 
Identification Program (UNALIP). Oye Ie all Vga owas 
intended to augment UNALIP Wiha ecm e get O8e 
identification routine using the work of Kogekar [1] but 
the existence of a well documented APL time series 
package [18] on the University of Alberta's MTS system 
made the inclusion of Kogekar's work redundant. 

The University of Alberta Identification Program is 
a modification of the Oxford Identification Package 


COALY )'e OXIP, developed for implementation on an IBM 
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1130 digital computer with 8K core memory capabilities, 
was obtained from Dr. D. W. Clarke of the Department of 
Engineering »Science,.University of Oxtond: «(Ox Go mds. 
England. UNALIP contains four identification algorithms 
along with various statistical and graphical options. 
The four identification methods are: 

1) Generalized Least Squares (STIME) 

2) Maximum Likelihood (ASTRM) 

3) Instrumental Variables (IVOFF) 

4) Unbiased Estimator (UBIAS) 

inesswUniversitye.sf ~-Alberta, Identification. Program 

attempts to identify the parameters of the pulse 
Deansner  s.unchion. lof. thersprocess. under. consideration. 
Tne. process, is.treated..as.a.'black box'. Therefore no 
DivswcalL ora mathematical ~description of .the system,is 
required, only the input/output data sequences are used. 
The source code is written in FORTRAN and the University 
of Alberta's Amdahl 470V/6 digital computer was used to 
implement the program. 

Chanter. ebasmanacutlune of (the etheory jofsthe four 
identification methods. The nomenclature and symbols of 
the methods were intentionally kept consistent amongst 
themselves to avoid confusion. 

A description of each of the major program modules 
was done in Chapter 3. An extensive overhaul of the 
program structure was required because the Amdahl 


computer's speed and size did not put limitations on the 
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3 
data handling of the program as was the case originally 
when OXIP was implemented on an IBM 1130 digital 
computer. Matrices and vectors were substituted for 
file handling, the maximum model order was increased, 
the mechanical switches which controlled program flow 
were replaced by software switches, improved program 
commenting was accomplished, simpler coding and a switch 
to FORTRAN IV were a few of the changes that were 
MistVigirceds resulting, in UNALIP. 

In Chapter 4, testing of the estimation routines 
Was carried out. The procedure used to determine the 
accuracy and speed of the various routines was the same 
as that used by Isermann [2]. Three different processes 
were tested using generated data of three different 
pseudo random binary sequence run lengths and two 
different noise/signal ratios. 

Chapter 5 contains the work done to estimate the 
pulse transfer functions for two pilot plant units. One 
UNnisAnewasmuacteaiquid-liquid, concentric tube heat 
exchanger. The other unit was a binary distillation 


column with eight bubble cap trays. 
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CHAPTER TWO 


THEORY OF IDENTIFICATION METHODS 


2.1 INTRODUCTION 
Time secolon describes= "tne “theorywror the-*four “off-line 
Pechulorecacion- alrgoritnms ~*used an?*the University of 
Alberta Identification Package. Most of the content of 
Pikomcnaptermirs, NOL: OT gifal Dut is*taken” from’ drfrrerent 
sources that deal with the techniques under 
consideration. The nomenclature and the symbols have 
been kept consistent for all techniques and special 
Sccencion has “been *"dinected "to clarifying certain 
aspects of the theory not considered in detail in the 


Oraginhal “sources, 
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2.2 STATEMENT OF THE PROBLEM 


Cousidér avidiscrete time, single input, single 
output system that can be described by the following 


linear, time invariant, difference equation. 
n 
Me Taha a) Cone) 


where: 1) fl is the noise free output of the system 


2) u, is the input of the system 

3) n is the order of the system 

4) k is the time delay of the system in terms of 
an integer number of sampling intervals 

53) ass Dd, are the system parameters 


Let the noise free output be disturbed by an independent 


zero mean additive noise signal. 
See ceev toe Gee 


where: 1) yy is the corrupted output of the system 


Za an is the independent zero mean additive noise 


signal 
The additive noise signal is considered to be the random 


process output from a general filter driven by a white 


mouse. source. This can be expressed as 
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where: 1) Ww, is the white noise source of the general 


Peicer 
2) nn is the order of the noise process (filter) 
3) Cis d. are the noise process parameters 


Equations (2.1) and (2.3) can be rewritten as 


[14A(z7')]v, ea ze eb zeeyus (2.4) 
[i*c(z ‘dle, = [1+D(z7') Jw, (2.5) 
where: 1) AC ze) - a,z +a at ee 
aera sper re biz fy yan Cae 
3) HAT OS = o427 + see + Cacia 
4) Dee) = d,z-' ae tera oe Ga wziaa 


Sleaze Se etnez—cranstorm evariable used! as; the 
backward shift operator (z7Ky, = Yt_i) 
-1 
Multiplying both sides of equation (2.2) by [1+A(z )] 


results in 


[14a(z7!)]y, = [1+A(z7')] vy, + bis(za, Ney €2 36) 


which can be rewritten as 


-1 -k -1 
[1+A(z ly, =+z¥"B(z du, + e, Ceca) 
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where e, = we Ut+aC27!)1ft4D(27 I /t14c(27') 1, 


Equation (2.7) can be rearranged as 


Y, = -A(z~ Weert we B( zi )u, +e (2.8) 

To «write equation (2.8) in standard vector matrix 
notation of regression analysis, consider a sequence of 
observations (u,y) to have been made of the system's 
input/output variables. Then assuming that bee for the 


system being considered 


W384, Be a, td Uy tee GU eee tea 
; (2.9) 
fee 84) i Git ot eon! tenon 
bu #...+D ou aaa 
1° t+N-1-k n~t+N-n-k t+N 


where N is the number of data pairs in the sequence 
Cu,y). Equation (2.9) can be written as 
k -1 


Va -ACz7 Dy +z Biz Ju+e (2.10) 


: = [ i 
Wheres 1) ¥ = LY, s---0¥i in 
af 
2) is CU, sees Up yy] 
iT 
3) e = Le, e+e sep yy 












es ar EP atewe < 
® a : ta > 
Re 
» 
+ ' ie y 
ms ' 


t's) 28ty, Be ane LLd Megs 


2 


(t.2° $+ use = + yuk 


+ 


; s a (eae 
kivzon tolaey faepoese 83 (6.99 ooldedos DESI BE 

/ (2 ‘ 
la 


‘Oo Semsupse 5 teptecoo.,ekeylene nejecemee 1, tee 
i ei ti 
s' weg A¥ Ro ce ovse iy. 


¥e 2 epee ‘ oF (¢ eu) see ihe 
oe) em ae 1 Oza j eee RTigaoe” = + fi T , oe ice IAS Vv. Pas 1 yi 
' + 
* : a : 
Lavopianos gephade 
ad ! 
o> m eB + ‘ees 4 Ud . ae. «3, Jee : 
4 i o a io ie n-a'n agt ee 
. ‘9 ¥ 
\ ~2 ' oe " 
~ hoy 44 af heed @f es 
_ ; t. Wie 9 " 1S tithe a 7 


SorHunse 94a ot e4ieq sIkb “to yd dpom. see. Be 
, é y 2% e4 7 
mm 


26 net@iaw 44 oso (0, 93- ‘csitsiga . 





Equation (2.10) can be expressed as 


= ite 


=] -kK-n 
Wier -z MV 4 Z Ufo, 2 uljal+ é 
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which can be written as 
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where: 1) 


ty 
[a,,---,a J 


A 
2% pb t= [b,,---,b] 
3) ae etn) Meee wean k 
X = 
RNs eas CaN ne 2 Sea gs PAV ERE oir irs pF 
a 
4) 6 = [a,b] 





2.3 ASSUMPTIONS 


In the statement of the problem, the following 

assumptions [3] were made: 

1) the functions [14A(z~')][14D(z7!)I/£14¢(27!)] 
and Cuenca Mi all have their zeros inside the 
unit circle 

2) there arel»no factors ‘common’ to’ all» three 


hain Benve? and bH.0 (254) ) 


polynomials [1+A(z_ 
3) for the type of system to be investigated, by =0 

and b,#0 
The assumption that the functions tes Wea. ei and 
Siege iilsD(zs)1/[14+C(20.)] have all their zeros 
imsiaqetethe ‘unit tcircle implies that the system is 
asymptotically stable. The second assumption implies 
that every state of the system is controllable either 


from u or Ww. For simplicity it is assumed that b)=0 and 


b, #0 ite’ 
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2.4 GENERALIZED LEAST SQUARES 
Generalized least squares is a technique of 
parameter estimation proposed by Clarke [7]. This 
technique attempts to produce unbiased estimates from 
Dowsye GCGata using least squares by filtering the 
input/output sequences until the bias is removed. on 
this section the least Squares technique is introduced 
and then Clarke's modifications are described. 
Thieveleastersquares estimation [5] problem dis the 


wove OumOi alesestimate O of 0 such that the loss 


function defined as 
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is minimized. 


where: 1) @ is the estimate of the parameter vector @ 


a ae 


ued 


is a weighting matrix, an indication of 
the degree of confidence placed in the 
individual measurements [29] 

3) e is the generalized error for the system 
under consideration and can be obtained from 
(ae 129 

4) V(@) is the loss function of classical least 


squares 
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The loss function is minimized by SCanmying out ‘a pariaal 


differentiation of V(@). 


“a ae XO) = 0 (233 5.) 
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8) | 
fees | = x 
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For the case of T taken as the identity matrix [6], then 


the least squares estimate of 9 is 
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Mie tyangespoune sides of ‘equation (2.12) by cx? xq7!x 
suoeusineg equation (2.13¢c) results in 
o 2 oteepe rt ys (2.14) 


Classical regression theory assumes that the 
elements of X are independent of the elements of e [7]. 
Therefore, due to the formulation of the problem, biased 
estimates will result because the elements of Xx are 
dependent on the elements of e (cf. equations (2.2), 
See weeand Ce st). s In this instance, the bias remains 
regardless of the number of data pairs [7]. This can be 
shown by rearranging equation (2.14) and taking the 


expectation of both sides of the resulting equation 
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From equation (2.7) and using a vector representation as 


in equation (2.10), e can be expressed as 


e = [1+A(z7')][14D(z7 |) Iw 
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Equation (2.16) can be rewritten as 
e = [14F(z7!)1w CAD 
where: 1) eS = jaa + ye +r 
2) f. are the coefficients of an infinite 


sequence produced by the PON ged yas Ot 
expansion of the expression 
[1+a(z~')]f14D(z~')]/014C(27')] 

Mnvseallows the iright hand side of equation (2.15) to be 


expressed as [7]: 
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Vico at Logeunecmnnners product of) .a.pair:».of 4s 


dimensional vectors 1 and m with the inner. product 
s 
Deeneovas: Sl ,m? Baiwed) ' lm. . 
Substituting the vector form of equation (2.7) into 


Coemecieit  wnand sesade, of equation (2.18) allows the 


expectation of the inner product of y and w to be 


expressed as 


E(<y weer el e/ze*atzm yt sey & ,W 
DiehCas oT ETA Cz) (2.19a) 


aqandstrom equation (2.16) it follows that 
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if it can be assumed that u and w are uncorrelated then 
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The right hand side of equation (2.21) can be expressed 
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* 2 * 2 
- fif,0 - f,f,590 e333 
nae * 2 
- f fio f CG —-.. 
E(xle) ign n+1° 1 n+2° 2 (2.24) 
~ OR ee ce 





The right hand side of equation (2.15) is non-zero 
unless the coefficients ek AY Gdn talt iee which occurs. Lit 
only when eave Wied (25). = Ree Sa which is not 
generally true. The estimates, 8, remain asymtotically 
Prased unless @ iS:an uncorrelated sequence. Note that 


even if in equation (2.2) is a white noise sequence, 


at 
the estimates of 9 are still biased because ey is not a 


white noise sequence [7]. 

If the coefficients of A(z~'), C(z~') and D(z‘) 
were known, then it would be possible to modify y and u 
Pyeraemoving average filter with coefficients of the 
autoregressive model of the noise process to give 


unbiased estimates of 9 [8]. The filtering process can 


be represented as 


hia Rotze lu (2.25) 
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vate sa Rotz))1y (2.26) 
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* 
where: 1) u is the filtered input of the system 
* 
Say pS ttem or itered corrupted output of- ‘the 
system 
1 -1 -IRO 


t= ro,z teeet TOR OZ 


4) ro, are the coefficients of the autoregressive 


3) Ro(z 


model of the noise process 
5) IRO is the order of the autoregressive model 
of the noise. process 
Peer eROTAS infinite’ ,i- thet’ coefficients tof s thet moving 


average filter are given by [8]: 
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(2.27) the following equation can be obtained [7]: 
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Since w is uncorrelated white noise, a least squares 
estimate using the filtered input/output is unbiased. 


The new least squares is 


Cons Oe 


* 


where: 1) is the filtered result of the matrix X 


X 
2) 3” is, the estimate of the parameter vector 9 
for the, filtered data 

Now the problem is one of finding a suitable moving 
average filter so that the input/output data can be 
modified to allow unbiased estimates to be obtained by a 
least squares estimator from this modified input/output 
data. This problem would be trivial if the coefficients 
of RoG2 ch) were available but the polynomial 


ee Cozen) and ner are not known. 


coefficients of A(z 
Clarke “(71 solves: this problem by designing a filter 
recursively until the final estimates are unbiased. ihe 


is assumed [8] that the estimates of the model error are 


a good representation of the noise sequence e. 
= - y (2-3) 
§ we yy, 3 


where: 1) sy is the estimate. of the model error 


(residuals) 
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2) a is the output from the calculated model of 
the estimation routine 
Using equation (2.27) and the definition of @, W Gan De 


expressed as 
w= [14Ro(z7')le (aon 
then using the assumption en $8, equation (2.32) becomes 
wo [14Ro(z27!)]g (2033) 


which can be expressed as 


§ = Gh + w (2.34) 
witere: 1) G-= ee eer 
2) n= lig 
h = Lro,,---r0r,5 
The least squares estimate of h is 
2 = eee (onBsD 


where h is the estimate of the parameter vector h. This 


estimate is unbiased because the elements of G are 


independent of w. The estimation of the autoregressive 
noise process model parameters allows the filtering of 


the input/output data to be performed. ticl@h bE Paaleat ego: 


nw “~ 


limit as h approaches h, the estimates 6 will be 
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20 
unbiased [9]. Since a finite approximation is used for 
the autoregressive noise process model parameters and § 
Was used for e, 9 is still biased and iterations must 
be employed to remove the bias. For each iteration, new 


* A % 
Peseoua ls 5 arelvevaluated using u,, y and 9 from the 


previous iteration using the relation 
BD NSy Sete ue (25530) 


where ae is the new estimate of the model error. 

The generalized least squares procedure is a 
"quasilinear formulation of a non-linear estimation 
problem" [8]. Convergence is not guaranteed unless 
restrictive assumptions are made about the closeness of 


the original parameter guesses to the actual values [8]. 
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2.5 MAXIMUM LIKELIHOOD 
The maximum likelihood estimation technique can be 
derived from the Bayes' estimator [9]. There are three 
assumptions used to derive the Bayes' estimator: 
a)) ther priobablityisdensiity fumetions: dof | cthe 
parameters ass dD. and f. are known 
b) the loss function which determines the cost of 
choosing the estimates of ass b, and fs is 
known 
ec) the probability density function of the noise 
sequence is known 
The maximum likelihood estimator drops assumptions (a) 
and (b), thus the parameters ass b. and f. are no longer 
random variables but known constant parameters. Astrom 
et al [3] have used the maximum likelihood estimator to 
identify» proces’) » parameterss/of»+systems*ecthat ‘can » be 
characterized by linear, time invariant, difference 
equations. ThiS, section deals with the maximum 
likelihood technique as employed by Astrom et al [4]. 
The linear, time invariant, difference equation can 
be obtained by using the vector forms of equations 
G2445 Mane? 27s) trand u€2:tiion atio rgivie 
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BBuze yD (2.37) 
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will be random variables that have a joint conditional 


probability density function represented as 


ply|u,a,b,f] (ase) 
PeesomIyephy|Urano stl tis thew joint conditional 
DVOUt Ut VeriensSivye Tunce lon, Of the soutput 
given the input, the system parameters and the 

noise process parameters 
2) SS eae 
A posteriori, the observed values of the random 


variable are known from direct measurement. 


a ! as ! 
Jom tothe taN 7! tan 


eremproblemeicseuhus sone.of determining a,b, and f£ from 


Vite s 

where: 1) a peuthemestimate sof the parameter vector a 
ied) b is the estimate of the parameter vector b 
33) f is the estimate of the parameter vector f 
0) $e ae Lyis tae wy) 


5 Jay, is the observed corrupted output of the 
system. 
The functional relationship between y' and er b and f is 
the same as the functional relationship between y and a, 
b and f [3]. This relationship is written as 


“a “aw 


Lly*|uya,d,f] (2.39) 
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Piers  Coalkled thew likelihood funetion. Usually 
Mattie |u.a,b,£)) is used [4] because the logarithmic 
function is monotonic and the maxima and minima of 
Miveru.a,b,f] occur at the same values as the maxima and 
fminama of plylu,a,b,f 1]. The logarithm of the likelihood 
function has the useful property [3] of continuous 
partial derivatives for all parameters. Since the 
Minima are finite, the gradient of the function will 
thus vanish at the minima. Therefore values of A b and 


aw 


f can be obtained by solving 


yew inliet iu, 2) oe) (2.40) 
3@ lg = @ 
T 
moere: 3) .0 ="la,0,.1°) 
2) g = La,b,f]2 


The solution must also satisfy 


a*. inbly'|u,@l| < 0 (2.41) 


90° Ig = @ 


Based on the assumption that the noise sequence w has a 
Gaussian distribution [4], the likelihood function of w 
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L{w] = 1 exp -W W ean yeh 
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Where o is the standard deviation of the noise Sequence. 
menwaneg the logarttnm of both sides’ of equation (2.42) 
results in 


In L[w] = -Nln2n - Nino - wi 


|= 
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The joint probability density function of YY sandeey "46 


dgerined as [10] 


pCy4 yo" = aioe lal PR (2.44) 
[ay| 
where: 1) y = [u,a,b,fl 


2) |aw/ay| is the absolute value of the Jacobian 
Of. ay 
pance y is a one to one transformation of w and the 


determinant of the identity matrix is unity, equation 


(2.44) becomes 
p(y,y) = ply.w) (2.45) 


and combined with the re fen ition of conadvutitoned 


probability [11], the following equation results: 
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p(y) 


goth: (2.46) 


p(y) 


If y and w are independent then [11]: 
ply,w) = ply) p(w) (2.47) 


Pumeisingeregcuattons (2.46) “and (2.47) gives the 


following expression 


Diyviay eae Gy) DCW) = pGw) (2.48) 


ply) 


NOWMEGeCUrIIng etOmeneay Likelihood function; it follows 


that 
InL[y'|y] = 1nLIw] 


= -Nln2r —- Nino oe ee (2.49) 
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The likelihood function is tnuss 2 function of the 
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parameters a, D, f. and Aa oac We ies) (ONC. all ol tea. 


conditions for equation (2.36). Ror ean Lio Cv. it 28 


assumed that the initial conditions are zero [4]. 
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As in generalized least squares, a loss function is 
introduced which gives a measure of parameter estimation 
accuracy. Theses loss «functionsrthat is. usediyean ebe 


expressed as 
V(@) = (ww)/2 (24.504 


ite wenset bea seen.» that.iminimizing e-the..losssafiunction is 
equivalent to maximizing iIn([y|y]) because equation 
(2.49) has its maximum value when (wow) /2 is a minimum 
Srnceeeconec. Other” two terms of the right hand side of 
€duation -(2.49) are constant and do not affect the 
Maramazarionuof in(ly|y):. The variance, bee, of the 
Motes sequence can be determined when G is known such 
enat VB) is a minimum by using the relationship 


ao Ns min(w w)/N CZ ee) 


where min(wiw) is the minimum value Of wiw. 


Tieepos! 1s to” Obtain a, b and f Such that V(@y is’ a 
minimum. This could be done analytically if the loss 


function had a simple dependence on a, b and f but 


examination of equations (2.36) and (2.50) leads to the 


|o 


conclusion that V(@) is a quadratic function of a and 
while’ the’ dependence’ of * f “is more complex [4]. 
Therefore an analytical solution cannot be obtained and 


an iterative approach must be sought. Due to the 
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Structure of the problem, a combined gradient Gauss- 
Newton technique where fast convergence is obtained 
through the use of second partial derivatives was 
originally employed by Astrom et al [3]. However, the 
method used in UNALIP is the modified Gauss-Newton- 
Marquardt algorithm [9] using a modified Levenburg 
parameter [12]. An update of the parameter vector using 
this method is calculated by the following equation 

soa 


t+1 iG 1 


) Cen cn) 


IQ> 
" 
|Q> 


Ag ty Oe Stee 
- (¥8g5 (9) 


where: 1) ¥ is the modified Levenburg parameter 


2) is the matrix of second partial 
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derivatives of the loss function 


32) S65 is the gradient of the loss function 


4) afk is the value of the estimated parameters at 
thesa't"thiteratron 
Since w is not available, an approximation is made 
to determine the noise source using the following 
relationship obtained from equation (2. 36) 
RPE eUEIeACz@ Sf = zece(z@ Jul/[i4F(zs)1 * ¢2.53) 


where: 1) w is the approximated noise source Of aatine 


model 


2) Wider Cosy is an autoregressive filter 
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Thus the loss function becomes 


(2.54) 


=> 
ct 


V(Q) = wow/2 =) / 


nm & 


1 


where We Leman, element of the noise) sequence w and the 
maximum likelihood estimate of the variance of the noise 
sequence becomes 


o° = min(wiw)/N (2.55) 


where ae is the estimate of the variance of the noise 
sequence. 
The partial derivatives of the loss function are 


obtained by differentiating equation (2.54) to give 
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where ®. is the "i"th element of the parameter vector Q. 
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2.6 INSTRUMENTAL VARIABLES 

Atssimpiie wolgtion (to).the Upnoblem sofiibiaslein 
parameter estimates due to observations contaminated 
with noise is the instrumental variable technique. The 
instrumental variable method is the same as the least 
squares technique except that the bias is eliminated by 
filtering equation (2.12) using the instrumental 
variable matrix. THps seetiometdeals with the 
instrumental variable technique as it is used in UNALIP. 
Using equation (2.12), the generalized error for the 


system under consideration can be written as 
e=y - XO (2.64) 


Unbiased estimates of 9 can be obtained from a least 
Squares solution if both sides of equation (2.64) are 


premultiplied by wi eiscre 


We = wiy ees (2. 059 


= 
WP< 
Ip 


where W is the instrumental variable matrix Which must 
satisfy the following conditions: 
a) ELW e] = 0 


b) W X is nonsingular 


If the elements of the instrumental variable matrix are 
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chosen to be uncorrelated with the residuals then 
ECW e]=0 and the problem reduces to a least Squares 
amalbysis with no bias [14]. Using a similar logs 


function as that described in section 2.4 


V(8) = (We) We/2 (2.66a) 


and applying the least squares method to this loss 


function results in the following equations 


SU een yt yen rake 0 
38 [0-6 

Wy = W'XO 

o = twexi7'w" (2. 66b) 


Pep CODLen Gi SerencCoun ered © in trying: toe find the 
elements of the instrumental variable matrix that will 
satisfy the stated conditions [14]. It has been shown 
Paap there, exists optimal instrumental Variables [15]. 
in UNALIP, the input sequence, u, and the calculated 
undisturbed output sequence y are used to construct the 
instrumental variable matrix. The calculated 


undisturbed output sequence is obtained by taking the 
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to be the parameters of the 
The instrumental variable matrix 
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The problem of bias in parameter estimates due to 
observations contaminated with noise is dealt with by 
introducing a new loss function to eliminate the bias. 
James et al [16] first show why the conventional least 
squares loss function can not handle the type of system 
being dealt with and then set forth a loss function 
which will handle the system. 

Using the generalized error, equation (2.64), the 
loss function can be defined as [16]: 

rT, 


V(6) = e ei= Cy — XO) 


Bape ay) (2.68) 


Conventional linear least squares estimates are obtained 


By minimizationfoft this loss function. 


0 (2.69) 


The authors [16] of this unbiased estimator method 
manipulate the loss function to show that the classic 
loss function is of no use for the type of system 
considered. This is done by taking the expectation of 


both sides of equation (2.68) to give 
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E(V(@)) = E((y - X@)"(y - x@)) (2.70) 


Equation (2.70) can be broken up into simpler parts by 


using equation (2.2) and the definition of X 


GCG) ERE G(ver men ECP 415 )G) (Vase (2211p 
aeeL eater ro )9.),) 
where: 1) 
oer) oe Ln Vie 3) ig Aer eee 
Bes : 
NN Wai ae Ne te Neeka sateen ene k 
2) 
mar ei, ats ace Lo eae) 
3S = 6 
SELON = Wee tan ony og ce 


Equation (2.71) can be manipulated in the following way 
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= E((v - Po)? (v - Pa)) + 2E((v - Pb) Ce 
- $0)) + E((e - $8)" Ce - $0)) 

= (v - PO)" (v - PG) + 2[E(v ec) - 87E(S'v) 
SERRA Haase Ses 
="207E(S7e)"+ BTETS's)0 (2.72) 


and remembering that the expectation of two uncorrelated 


sequences is zero [17], the following equations can be 


written 
E(v'e) = 0 
E(P'v) = 0 
E(Pte) = 0 
E(P'§) = 0 


where: 1) O is the null vector 
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So equation (2.72) can be rewritten as 


E(V(@)) = (v - Po 


=Oge Sie) + Een €) ACESS, 


Micmtirsoe cverm “of equation (2./3)) is the conventional 
linear least squares loss function used when noise free 
data is available. The second and third terms lead to 
the bias in the estimates while the last term is related 
to the variance of the estimates [16]. 

The method of James et al [16] is to reject the 
conventional loss function, equation (2.68), and use the 
following loss funetion to eliminate the effects of the 


Second and third terms of equation (2.73). 


J(9) = (y - X98) Cy - X98) - 9 RO + 28 Q C2et a) 


Mere: 1) R = 
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The estimate of the parameter vector 9 is then 


determined by classical least squares [5] to give: 
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CHAPTER THREE 


ORGANIZATION AND STRUCTURE OF THE UNIVERSITY OF ALBERTA 


IDENTIFICATION PACKAGE 


3.1 INTRODUCTION 


This chapter details the purpose, method, flow 
diagrams, box diagrams, variable descriptions and 
Subroutines called by the major program modules of the 
Unmaversdity sot Atberta thdentificationashdackage’. In 
general, the package accepts input/output data from 
noisy, single input/output systems, and estimates the 
Mest discrete, pulse transfer function models for the 
system and the noise processes. It can handle systems 
with time delays, integrations and nonstationary noise 


processes. The model obtained is of the form 


y = a tates yu + [14D (27) Jw C31) 


Eira(zpey] Glican a. 
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where: 1) 


2) 


Sy) 


ni) 
INT = IYDIF, —“IUDIF 


IUDIF is the number of differencings of the 
input sequence 
IYDIF "Ys “the number of Gsiterencing’s *of= tne 


Output sequence 


The optimal model must pass the following tests: 


1) 


oy) 


3) 


4) 


ay) 


There 


variables 


mean square error improvement over a simpler 
model 
Significance of individual parameters 
nonsingularity of the parameter covariance 
matrix 
uncorrelated residual sequences 
zero cross correlations between the input 
Signal and the error signal 

are certain restrictions placed on various 


duew lon naty 1 set eVeCLOT, 8S176 5 and) elogical 


requirements of the FORTRAN language. The restrictions 


ares 


1) 


2) 


3) 


Che aAnOumtetO imo a Valet oo eres tracted: to (550 
input/output pairs 

the model order of the process must be between 
One. and tive torso liMe {3.9}, IVOFF {3.11}~and 
UBIAS {3.12} and must be between one and three 
for ASTRM (3.103 

the time delay must be greater than or equal 


to zero 


All programs are written in FORTRAN IV. There are 


seventeen major programs all interlinked by a steering 
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program, STEER {3.3}, which provides the user with the 
aoLlacy, FO. control program flow. Each time program 
control is transferred to a new major module, the 
program's name is written and underlined. Normally, 
program flow is automatic and the user does not need to 
take personal step by step control. A flow diagram of 


tuew "normal" mode of control is: 


DAGEN or DATGET 


mworoptions of note are. 


1) the MONTE CARLO facility which allows the user 
to run the same data with the same model order 
and time delay range through all the 
estimation routines without hav.ingveto 420 


through the entire eycle, theta eis. one 
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subroutinesmnbior {3.14}, °NOISE {3.15}, -FINAL 
i310) sanaerpior, {3.17} are skipped. 
the subroutine CRCOR {3.8} computes the cross 
correlation function between the modified 
input/output sequences. If the input sequence 
is a pseudo random binary sequence then an 
unbiased estimate of the weighting sequence is 
produced. [If the input sequence is not a 
Dseudo random binary sequence the cross 
POtmecaceole rune. On OWL leenave the” correct 
Shape but the amplitudes will be incorrect. 
TivseecogTeclatnhoneeanely Sis <i1S° ._provided “by 
UNAeireepae oCiis —faciiaty 2s not sexamined in 


Cha seones.1-s < 



















‘oo = ' - —— oe “| 
: ° The, ’ “we - dak ee ' a 
‘ ~ 7 
: , , 7 
fl ’ os 
ra , —§ io 
i. y - _ @ 


iP AA F > " > ae x," ; eae 
JavTH. CELE) RREOW, STEELE) Gel2e senna 


d a ry . r dh 
baactie era {Pe NeIG% ovatae 
- ; ¢ 
oo (8.8) MOPRD snrxuoraer ss 
i : io 
»sPlibom easy asewied anpatoaut no2iszerst 


- 7 
40s 
_ 


fw 


. © , 
faune2 Judol sag 7] .essifevpses rugsve\ewgak 
a 


fant 6oneupse yurBnid, Mothie: Gnmnseq. ¢ =} 


¢ ton -2% sensuoce. Cupdi gg 7 : spubp td 
22D st gane.p¢se veaed mohjat gbueag 
5 and SVvex ‘yw notdon? notgveleaides 

at’ 'itiw epoddifewms ia ‘ad ogee 

Pyanng J tien leant TolitsLe tos ewePrt 
‘titoe? etd? dua FTISAAY 


e, @hasnz eid 


ae ra 


42 


PROGRAM MODULES 


Bun Pur: 


the mainline program. OXIP initializes all 
major system variables and transfers control 
Sop rograme sil OW tO, Lhe subroutine To DEER: 
AS) s The flow diagram for the OXIP module 


is figure 3-1. 
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VARIABLES: 


COVM 


COVV 


el 


Be 


IFILE 


Can D real vectors of length five used to store. 


the coefficients of the A, By, G, and D 


polynomials. 


the covariance matrix used to store the matrix 


product Cem 


the covariance vector used to store the vector 


product xly 


a real vector used to store the noise added to 
the unmodified output data in the subroutine 


DAGEN {3.6} 


a imealdamrectorjisugediitoy stonesathe residuals 
calculated in the estimation routine 
emda lE testa —nek¥e whene: Y2seis,sanwdntermediate 


calculation of the system's output 


a real vector used to store the prediction errors 
(final residuals) calculated in the estimation 


routuneseies) B2eq Kh —atle 


an integer variable that contains the logical 
unit number of an output file used to _ store 


calculated results from many subroutines 
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an integer variable that contains the logical 
unit number of an output file used to store the 


BeateVeCcCORem heart, 2, UyeUt. U2, Y, Y1 and, Ye 


an integer variable that contains the logical 


unit number of the input device 


an integer variable that contains the logical 


unit number of the output device 


an integer variable used as the major steering 


parameter that dictates program flow in STEER 


{3.3} and other subroutines 


an integer variable used as a switch to indicate 


if the Monte Carlo facility is needed 


a real vector that stores the parameters of an 


autoregressive model of the residuals 


a real vector used to store unmodified input data 


a real vector used to store input data modified 


by the subroutine DMOVE {3.7} 


a real vector used to store modified input data 


that is manipulated by the estimation routines 
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vector used to store unmodified output 


vector used to store output data modified 


subroutine DMOVE {3.7} 


vector used to store) the prediction of the 


output by the estimation routines 
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SUBROUTINES CALLED: 
eens. 3S} )06U6valows “linkage to the rest of the 


subroutines 
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Bao SUBROUTINE: STEER 


PURPOSE: 


METHOD: 


a general steering program that regulates the 
PLOW, BOtemexecution of “UNAL IP. The flow 
diagrams for the STEER module are figures 3-2, 


3-3 and 3-4. 


the flow of execution is determined by the 
system variables ISWCH, JDAT, JSWCH and MONTE. 
These variables are set by different 
Subroutines or by the user. This subroutine 
also sets up Monte Carlo runs and allows the 


user to exit from the package. 
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Figure 3-2 
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VARIABLES: 


IRO 


JDAT 


JSWCH 


NL 


an integer variable that contains the current 


autoregressive model order of the residuals 


an integer variable used in computed GO TO 
statements to determine which estimation routine 


is to be used 


an integer variable used as a switch to indicate 
the direction of program flow determined in the 


minor subroutines DGSW and DGSW1 


an integer variable used as a switch in: 
1) the minor subroutine PARSP to determine the 
type of calculation to be done 
2) the estimation routines to determine if 
ir Sc guesses of the polynomial 
coefficients and other variables were used 


in the subroutine STEER {3.3} 


ayeminceger® Variable. used  to® hold® the current 


system model order 


integer variable used to hold the lower bound of 


the system model orders to be investigated 
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integer variable used to hold the upper bound of 


the system model orders to be investigated 


an integer variable used to indicate program flow 
determined by the minor subroutines DGSW = and 


DGSW 1 
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SUBROUTINES CALLED: 


ASTRM 


CRCOR 


DAGEN 


DATGET 


DGSW 


DIAGN 


DMO VE 


FDISP 


{3.10} an estimation routine using the maximum 


likelihood algorithm of Astrom et al [3] 


{3.8} computes the cross correlation function 
between the modified input and the modified 


output data 


{3.6} generates output data for a noisy, discrete 
time system excited by a predetermined input 


sequence 


{3.5} transfers input/output data from user files 


to system vectors 


Supplies messages and program flow control at 


decision points in the program 
{3.13} takes the current calculated system and 
performs statistical tests to determine the 


"goodness of fit' with the real system 


{3.7} modifies raw input/output data and certain 


system variables are also obtained here 


{3.17} supplies graphical displays of data files 
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{3.16} gives various tabulated results for the 


system considered 


allows the.user to use free format input and 


provides certain error checks for the input [21] 


{3.4} describes various input variables required 


by the major program modules 


isneisseean Sestimation routine wsang: san) off—line 


(non-recursive) instrumental variable method 


{3.14} allows the user to determine the impulse 


response of the optimum system model 


{3.15} estimates noise models from uncorrelated 
measurement error data from the subroutine MDISP 


Ste as 


writes out repetitive headings 


{3.9} an estimation routine using the generalized 


least squares algorithm of Clarke [7] 


{3.12} an estimation routine using an unbiased 


algorithm [16] based on Slutsky's Theorem 
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3.4 SUBROUTINE: HELP 


PUT rOSE: to describe the various input variables 


required by the major program modules. The 
fpows diagram for sche HELP module is figure 


3-9. 
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ABLES: 


an integer variable used to determine program 


flow § as reortat 
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SUBROUTINES CALLED: 
FREAD allows the user to use free format input and 


provides certain error checks for the input [21] 
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3.5 SUBROUTINE: DATGET 


EURPOSE: 


METHOD: 


to transfer input/output data from user files 
to system vectors. The flow diagram for the 


DAGET module is figure 3-6. 


the input/output data is stored in two files 
DATIN and DATOUT that are assigned the FORTRAN 
logical unit numbers 18 and 19. The minor 
subroutine FREAD [21] is used to transfer the 
data from the files to the system vectors U 
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VARIABLES: 


ICARD 


INDAT, 


JJ 


LENGTH 


NUMBER 


USAVE, 


an integer variable that holds the number of 


lines in the input/output data files 


IOUTDA integer variables that contain the 
logical unit numbers of the input/output data 


sources 


an integer variable that is used to transfer the 


Gdavansirom the files to’ the vectors 


an integer variable that contains the total 
number of input/output pairs transferred from the 


files 


an integer variable that holds the number 
Of avalues On slay lane of the input/output data 


files 


YSAVE real vectors used to store one line 


of the input/output data files 
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SUBROUTINES CALLED: 
Papa “allows the user to use free format input and 


provides certain error checks for the input [21] 
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3.6 SUBROUTINE: DAGEN 


PURPOSE; 


METHOD: 


y= 


MoAv (za Il oz 


where: 


1) 
2) 
Sy 
4 ) 
aa 
6) 


te) 


8 ) 


to generate output data for a noisy, discrete 


time system excited by a predetermined input 
sequence. The flow diagrams for the DAGEN 
module are figures 3-8 and 3-9 and the box 


diagram for the process is figure 3-7. 


the output data is generated by the model 
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Ug Sa. DreceruerminedmetnpUL sequence 


z is the z-transform variable 
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w is a normally esd2tstributed) uncorrelated 


random sequence (white noise) 


y is the resultant output data sequence 
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baer aids, aa and C.de are the system 
i i i i 
parameters supplied by the user 
IDRS is the number of integrations to be 
performed on the output sequence 


IDRN is the number of integrations to be 


performed on the noise sequence 
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VARIABLES: 


AD, BD, 


EMS 


IDRN, 


IJ 


ZOUT.P 


IX 


JFRST 


JLAST 


CDy¥ aDD real vectors of length five used to 
svore=the coefficients of the AD, BD, CD esanruseuD 


polynomials 


a real variable that contains the mean square 


value of the generated noise sequence 


IDRS integer variables that contain the number 


of integrations to be performed on the noise and 


output sequences 


an integer variable used in DO loops to compute 


the output and noise sequences 


an integer vector that holds the output of the 


pseudo random binary sequence generator 


an non-zero odd integer variable that is the 
SY ic daw seed " of ** the** Gaussian and Suniform 


distribution generators [19] 


an integer variable that contains the number of 


the first input/output pair considered 


an integer variable that contains the number of 


the last input/output pair considered 
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1) an integer variable usedinasthan # cinirteian 
value in a DO loop to calculate 
y=BD(z~!)u/[14AD(z7')] 
2) an integer variable usiedijq@asthan ¢ eindétidial 
value in a DO loop Ke calculate 


e=[14DD(2~!)w/L14CD(27! 


al 
an integer variable that contains the number of 
input/output pairs generated plus the number of 


time delays plus one 


an integer variable that is initially equated to 
PA SOM that IX is. not destroyed by the minor 


subroutines GAUSS and RANDU 


an integer variable that contains the number of 


time delays in the generated system 


an integer variable used in the minor subroutine 
RANDU to return. the next value of the seed for 


the generator 


an integer variable that determines the length 
of the psuedo random binary sequence 


ie.) LENGTH=2 **NCODE-1 
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an integer variable that contains the generated 


system model order 


an integer variable that contains the generated 


noise model order 


an integer variable that contains the number of 


input/output pairs generated 


a real variable that contains the period (usually 
ijn tiples of the sample =time) of the 


Sinusoidal disturbance that generates U 


ae realise variable wused to store the white noise 


sequence 


a real variable used to store intermediate 


Se Oa PL Ol se.0 tan an Gi 


Area) | Varin Deemeoulaw ‘COMmvAins. the ) svandard 
deviation of the Gaussian distribution generated 


to give w 


YX DC real variables used to add constant values 


to all elements of the input and output sequences 


generated 
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YMS a real variable that contains the mean square 


value of the generated output sequence 


ZNSR a real variable that contains the noise/signal 


root mean square ratio 
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SUBROUTINES CALLED: 


DGSW 


FREAD 


GAUSS 


PRBS 


RANDU 


Supplies messages and program flow control at 


decision points in the program 


alliows ether user sitoeuse -fneereformat sinputr. and 


provides certain error checks for the input [21] 


generates a normally distributed uncorrelated 


random sequence (white noise) [19] 


generates a pseudo random binary sequence of 


variable length 


generates a uniformly distributed random sequence 
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PURPOSE: to enter parameters required by estimation 


routines and to perform transformation of the 


raw data to mean free modified data. The flow 


diagram for the DMOVE module is figure 


3-11 


and the box diagram for the process is figure 


3-10. 


METHOD: the data is transformed by differencing 


high pass filters 
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VARIABLES: 


AU, AY real variables that contain the means of 


the raw input and output sequences 


mu. AYY real variables that contain the means of 


the modified input and output sequences 


CRITE a real variable that contains the convergence 


criterion used in the estimation routines 


PORU, IDRY integer variables that contain the order of 
theitmighsrpass ij dfiliterssefor.othe Gnput ltandroutput 


sequences 


IFRST an integer variable that contains the number of 
Cheviiirsteetinput7outpute® paireitor*be used by the 


estimation and other routines 


Pow, Lecy integer variables used to calculate the 
camenecons tantstis=-spacedesafor the high pass 


filters of the input and output sequences 


ITMAX an integer variable that contains the maximum 


numberseof iterations used in the estimation 


routines 


IUDIF, IYDIF integer variables that contain the number 
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78 
of differencings of the impute! ‘and output 


sequences 


JFRST1 an integer variable that contains the number of 
the first input/output pair to be modified by the 


Subroutine DMOVE {3.7} 


KFRST an integer variable that is used to initalize the 
first values of a data sequence that is to be 


operated on by a low pass filter 


KL an integer variable used to hold the lower bound 


of the time delays to be considered 


KU an integer variable used to hold the upper bound 


of the time delays to be considered 


we Wo an integer variable that is used to hold the 
number of the last data pair to be used by the 


estimation and other routines 


RMSU, RMSY real variables that contain the root mean 


squared values of the modified input and output 


sequences 


mous TCY real variables @thates.comeain “thes wine 


constants for the high pass filters of the input 
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where: 1) T is the sampling interval (equated to 
unity for this package) 
even cCeuls the time constant of ‘the filter 


Gaquatedsto LTC) 


UMEAN, YMEAN real variables that contain the sum of the 
means of the raw and modified input and output 


sequences 
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SUBROUTINES CALLED: 
FREAD allows the user to use free format input and 


provides certain error checks for the input [21] 
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3.8 SUBROUTINE: CRCOR 


PURPOSE: 


METHOD: 


to compute the cross correlation function 
between the modified input and modified 
output. The flow diagram for the CRCOR module 


is figure 3-12. 


if the modified input sequence is a pseudo 
random binary sequence, the scaling is chosen 
to result in an unbiased estimate of the 
system's weighting sequence by the cross 
correlatvon=tunction. If the Modimied input 
sequence iS not a psuedo random binary 
sequence, the cross correlation function will 
have the correct shape of the weighting 
Sequence but’ the amplitudes will be wrong. 
The mean of the weighting sequence is removed 
from the sequence and a graph of the cross 


correlation function can be obtained. 
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VARIABLES: 


INOC 


JLAST 


JTOTC 


LEN 


NCODE 


NCORR 


an integer variable that contains the number of 
times the pseudo random binary sequence is 


repeated 


an integer variable that contains the number of 
the last input/output pair considered Dye oat his 


program 


ateeenveger variable “that’ contains the total 
number of input/output pairs considered by this 


program 


an integer variable that contains the length of 


the pseudo random binary sequence 


an integer variable that determines the length 
of the pseudo random binary sequence 


ie.) LEN=2**NCODE-1 


an integer variable that contains the number of 


eross correlations to be performed 


1) a real variable used to determine the two 


different levels in the modified pseudo random 


binary sequence 


is 
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2)" a real variable used to determine the mean of 


the cross correlation function 


a real variable that contains the scaling factor 


forl.cne crossecorrelation function 


a real variable that contains the difference in 
the two levels of the modified pseudo random 


binary sequence 


a real vector that holds the cross correlation 


function 


real vectors that contain the input and 
output sequences used to generate the cross 


correlation function 


a real vector and array that are used to 


pass functions to the minor subroutine GRAPH 
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SUBROUTINES CALLED: 


DGSW 


FREAD 


GRAPH 


Supplies messages and program flow control at 


decision points in the program 


allows the user to use free format input and 


provides certain error checks for the input [21] 


allows the user to plot on a hard copy terminal, 


a graph" ofSa function supplied! by -the “user 
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STIME 


— 


PUR POSE: 


METHOD: 


to calculate the coefficients of a model of 


given order, choosing the best time delay over 


the allowed range, using the generalized least 


Squares method of Clarke [7 1. The flow 
diagrams for the STIME module are figures 
3-13, 3-14 and 3-15. 

unbiased parameter estimates of a noisy 


process using a least squares solution can be 


obtained by Picure. in the Inputs OUbDUL 


sequences using a moving average filter. The 


coefficients of the moving average filter are 


the coefficients of the autoregressive model 
of the noise process. This, leads: Stomethe 
estimate 
a % *T *_-1_*T * 
oiee= 1 ot X Bi } We oe Sipe) 
® 
where: 1) X iSurtienenincered result /ofeune 


matrix X 
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% 
® is the estimate of the parameter 


vector @ for the filtered data 
% 
Yo ols the filtered corrupted output 


of the system 


ie Sy a | Wak: | 







e 
7 ay Pps 
\ofemn1eg sty Jo s2enbeee eae si # 


, . ood : hy ~*~ 2 _ 
adeb bavec fs? od, WP @ v9too" j 
: Pe : ‘> 
‘ugiwo boettuston bets t 289 Bay wa ce 


ngfege sd? 14 


~ 





88 


ad 
a 

oe 
v 


Peers pgrcocccoces 


soe 


eerse eevee eoeee eoove @vece eeose 





y E F 33 
a5 | ice: eee SPE Sab UB eae eee AE 3 | Bal 
x | EER F bE Fea 
a4 a an : 3 4 BH ral 
a eee aS ey hd fer Ase PCR 
He 
a ens hh. ois Sie eae 
E 
es F 
8 eee ee ee @eeoeve 
@ z as mas i re | : g | wa eee 
coe as a d ap d eal 6 cab FE 3 Be chats pete coe eS BL eee Hg ced 5 | 
3 ae ae b, f. a, R F* : war 
RE Ra) Bs oe: oa Lea Be 


ee ie ieee ies 
le Te EB ee poe. 


Flow Diagram for STIME 


3-13 


Figure 


ar) 
oy om, 
ie 
. 
er gente . 
‘ } . re, 
- ears vi we OS, 
Se eer, 2 * eee eee: 
wads o- mie i 
¥ 
ng ee er nee 
irony <3 . ne 
a ot sealed TS ae’ 
= bee? eee b Vanes 
- * 
a 
‘ 
nd 
—a— = = £ 
; a 4 a* s 
vr ot O18 7 Ow: 
= AA e 
=: r 
_ na — 
a ee weer 


4 
eT, 

















° 
a - = an, — 
. . . * + 4 * 
°* > . * alee: ai Z 
‘ aie 44 - - ve * « . . 
* ‘ * . - 
axmrmnioe — — ave 
a} ive 
: : rf 
ae 7, o 
« al 
* MT Dr A 








hd 

M 
ae ey 
. ’ ° ° 
ue a = . aRITE . 
©. FErmrsee? = *————— > paaaerTERS ¢ 
.. . . . 
e a) . . 
ie at pereeeecceceenees 

oe.) 
< 
v 





. 
et 
owes 
ae. e 
° oe age 
. 
coos 
pra lita “arate 
e 
e 
° 
Py 
ececcccrecccccese 
e 
Seanad 4 beneeemeees 
° ° 
° sare CouexT * 
o vaniaBlLES ° 
° Py 
- ° 
seuccccsaccccsccs 
v 
o% 
Pi eo. see f 2a 
g@ a a . 
e . wT . . 
© LPDIFPOCRITE? -¢—————»* LCRENEST TEES? 
ie: .* . . 
:. 4) eevee cre cocceree 
a) a 
’ v 
e ° ° ° 
° mpc : ° 
° CONVERGENCE * ° so noe ° 
. ° ° ° 
Ps e . e 
ccocccesecccoccse oeececsccecesscce 
sone : 
. 2 | 
ons : 
. i “ 
ooo : v 
° 





. . 
‘ae 
. . 
toes 
’ 
. - 
oe “ 
° °. . . 
e * res . ‘or 
© SIGINCTRETAD | Cem 0 ti be 2? taee 
. . e 
“* ° . . 
. . 6 . 
a.) ° 00 cove 
o 
eeee 
v 
ao ‘ 
p -. a Ory seweenseccercoccs 
. ., ° . . - 
. * wo : ° © Inchras= cre * 
= liaTTei? ‘——— * EFESTT STEP SIZ = Slat = 
: . . . . . 
. . ° . . . 
ee eeeceeccecseeeeees peccecc ccc ccocccs 








sesie 
° pe 
eceee tree ececeene oe 
. e e 
© Ser inv a * 
° imc . 
. ° 
. e 
’ 
eores eee een ene 
. ° 
° . 
° UCKEXENT JTEES? 
. 
. . 
woceeccccccccces 
“e 
eeeeegpeceeceeses 
°. e . 
e ° . . 
°. JTEES!? —_—_—° WEITZ sea o e 
° 
ace . . . 
wreteceereccecces 
a "i 





Figure 3-14 


Flow Diagram 





for STIME 


89 


coorerscccowccons 


SET (1aTT 


eens 
eonee 


= 
_ 
= « 
or ~~ 
* 
+. ~ é 
——— 
- a 
*, — 
a Or 
— e--~ 
 t 
e, 2. 
- 
= « _ 
’ 
om ae 
- 
. 










® 
SERESTORE PSes lame 
© VARIABLES OF °¢ 
SS ConkeENt Ko ° 


sins 
"load e 
oe °. xO 
© THETAVTHTA? .° 
Ad ° 
oa a 
ee 
SeYES 
< 
v 
seeeor) cecccoooce 
= e 
e e 
=> (NCREMENT KE ° 
» ° 
ba e 


’ 


2 e 
© RESTORE 35ST ° 
e SYSTE4 s 
© PARAMETERS be 
e 


Figure 3-15 











2 ° 
© CALCTLATE ° 
© PREDICTION e 
e ERRORS oO 
se 


veteantetteeerees 
of e 
* WRITE “AJOR °¢ 
: VECTORS e 
< ° 
° e 


Se eee 


Flow Diagram for STIME 


90 


me. 
5 





91 
VARIABLES: 


Paaeeo ene UNnNeyDETGeEESHILRO; KKg RRO, SIGM, UU, YY are 
real sand integer variables, vectors and arrays 
that. store the best values of A, B, COVM, COVV, 
PEvoreeral,e re enO, K, RO, “SIGMA, U2 AND Ye 
respectively for the model under consideration by 


the estimation routine 


eeeeoi ee CKOVYM,. CKOVV, DDETE, ROR, SSIGMR are real 
variables, vectors and arrays that store the best 
uaLuesPos Ar, “Bit LGOVM ye COVVTIDETE.® RO ean'de SIGMA 


forseach iteration of the estimation routine 


ADIFF, BDIFF°real vectors that contain the difference 
between the previous iteration and least squares 


estimatve of*A and B 


ANEW, BNEW wealeevectons tthat scontain )ithesar parameter 


updates of A and B 


ASAVE, BSAVE real vectors that contain Aean dieb af rom “the 


previous iteration 


DETE a real variable that contains the determinant of 


the covariance matrix 
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a real variable that contains the step size for 


the approximation of A and B 


a real variable that contains the average of the 


residuals 


an integer variable that: indicates if the step 


size has been halved 


an integer variable that indicates the type of 


exit from an estimation routine for a particular 


time delay 
1 = convergence 
2 = divergence 


3 ITMAX has been exceeded 


an integer variable that indicates if the step 


size has been reset to an earlier better value 


an integer variable tnace contains) .the -current 
number of iterations in the major loop of the 


estimation routine 


an integer variable that indicates if the step 


size has been increased 
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an integer variable that contains the dimension 


of the square matrix COVM 


an integer variable that limits the number of 


times the step size can be halved 


an integer variable that is used to hold the 
value of the current time delay being considered 


by the estimation routine 


ate integer variable that indieates' the type of 
exit from an estimation routine for the optimal 
time delay 


1 


convergence 


2 = divergence 
3 = ITMAX has been exceeded 
area levariable that contains the convergence 


criterion calculated by the estimation routine 


real variables used for intermediate calculations 


of ANEW and BNEW 


a real variable used to contain the mean square 


value of the residual 
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THETA areal variable used to determine if the current 


iteration is better than the last iteration 


VTHTA a real variable used to determine which time 


delay gave the best results 









ae 7 ‘an i4 Pan 


’ - Re 


eins ey “As snkeceta te enw aie a4 e% 


y 
‘noeige ysis! seal: adi ans “ed od: an me vs 
ii n : 


ie fae WA Sem Sey a 


t f a 
. ; — a; 

: ' an | 
iw. solersveb of Seu olde hwy,  See- 
; - 7 
is) i, 


+ I1gee4 peed acz “vas 


" ' 
< . - . 
,; : an de ; ° , Pe Ee > ad Fr 
| y 3 +e, Lt ear) Sp , Mera A 
J i ree - A ; we +; 


SWITCH 


IREST 


IIATT 


IHALV 


oe 


POSITIONS: 


1 


re 


" 
= 


Wt 
ass 


initial position of the switch 
indicates that the step size has been 


returned to a value that gave better results 


initial position of the switch 
indicates... that.»-thes step size -chas, been 


increased 


PniGualtpositironnoef athe lswiteh 


indicates that the step size has been halved 


acy r* 


tad 


re il) 


- My are 
ssi2e gave od sada codon bin: 


(es a “ 
fod ovey todd) Sedav 2 of bonwene 
ie ne 
: * _— > i 
we.) 
inztiwa end tm aeivicog feittal f & 
j 7 J 


ssee ‘cede sh? dandy esteotbat 














| i -— 2 
ey, } yy 
we le 
+f Ta 
tos twa aad Me get) Laed fpidia: 


te 
Rar 
7 


¢ 
-_ 
) 


i 


xz oft Yo s@bddegg Isisiaz 05 


= rs 
>is qsate odd des’ eeteolpad S 





96 


SUBROUTINES CALLED: 


AFILT 


FCURG 


MCOMP 


PARSP 


PREDS 


PRINT 


RUNAV 


filters an array using an autoregressive filter 


computes the parameters of an autoregressive 


fixed order model 


computes the covariance matrix, the covariance 


vector and the determinant of the covariance 


matrix 


prints’ out the current. model order results 


computes the deterministic predictions of the 


discrete model 


prints out repetitive headings 


filters an array using a moving average filter 
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3.10 SUBROUTINE: ASTRM 


PURPOSE: 


METHOD: 


to calculate the coefficients of a model of 
given order, choosing the best time delay over 
the allowed range, using Astrom's et al 
maximum likelihood algorithm [3]. The flow 
diagrams for the ASTRM module are figures 


3-16, 3-17 and 3-18. 


£60 Maximize the likelihood function, the loss 
function 


“a 


1 N 
hae |W C3743) 
etal t 


V(G) = 


is minimized by use of the modified Gauss- 
Newton-Marquardt algorithm, using a modified 
Levenburg parameter. Thus the parameter 
update equation is 

“t+ et 


eh iys 7800 -)17083(0.0 (3.5) 


00. 
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So 


where: 1) @ is the vector of estimated 


parameters a, b and f 


2:) S 5 is the gradient of the loss 


function 
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S30 is the approximated matrix of 
second order partial derivatives of 
the loss function 
y is the modified Levenberg 
parameter which multiplies the 


diagonal elements of 33g by a 
constant 


W is the approximated noise source 


Ve 
of the model 
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VARIABLES: 


A2 


BESTT 


DEL 


DELI 


DELUPA 


eR P 


ee OL 


IUP 


a real vector used to contain the best residual 


sequence of a particular iteration 


a real variable used to determine if the current 


value of SIGMR is better than the previous value 


a real variable that contains the step size for 


a change in the parameters 


a real variable used to modify the step size DEL 


aarreal variablevy that recontains Ci Cmat Lise 


trial point for the step size of a new iteration 


an integer variable that contains the number of 
times the modified Levenburg parameter ok 
So epeved soy.) aneconstant § for. - ay particuder 


Drernation 


anmeunmtceger variable that contains the total 
number of times the modified Levenburg parameter 


is multiplied by a constant 


an integer variable that contains the total 


number of times the step size is successfully 


doubled 
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an integer variable that contains the dimension 
Gr’. the approximated second order partial 


derivative matrix 


an integer variable that contains the current 


dimension of the square matrix COVM 


a real variable containing the modified Levenburg 


parameter 


a real square matrix containing either the least 
Square matrix (x? x) or the approximated second 


order partial derivative matrix 


SWITCH 


IDOUB 


IITATT 
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POSITIONS: 
PObinatial.~position, of.the, switch 
yPeuUmeindicates, that the step size has been 


increased by a certain value 


oO 


invcialweposition of the switeh 
= 1 indicates that increasing the step size has 
given estimates that are poorer than the last 


estimates 
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SUBROUTINES CALLED: 


AFILT 


BCOMP 


CCOMP 


INVERT 


PARSP 


PREDS 


PRINT 


RESID 
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filters an array using an autoregressive filter 


computes the matrix of second order partial 


derivatives of the loss function 


computes the gradient of the loss function 


computes the inversion of a given matrix 


prints out the current model order results 


computes the deterministic predictions of 


discrete model 


prints out repetitive headings 


computes the residuals of the estimated model 


the 
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3.11 SUBROUTINE: “IVOFE 


PURPOSE: 


METHOD: 


tomcaliculate= the» coefficients of a model of 
given order, choosing the best time delay over 
the allowed range, using the instrumental 
variable technique. The flow diagrams for the 


IVOFF module are figures 3-19 and 3-20. 


unbiased estimates of a least squares solution 
ean be obtained by premultiplication of 
Srnivasuon G2. by Wo the instrumental 


variable matrix. This leads to the estimate 
x7 wy (3.6) 


where: 1) is the output sequence 


n= kK 


2) is the instrumental matrix 


is the input/output matrix 


3) 
4) 


|D > 11>< 


is the vector of estimated 


parameters 
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SUBROUTINES CALLED: 


BCOM P 


CCOMP 


PREDS 


PRINT 


RCURG 


RUNAV 


computes the instrumental matrix 


computes the instrumental vector 


computes the deterministic predictions. of the 


discrete model 


prints out repetitive headings 


computes the parameters of an autoregressive 


model by a recursive procedure 


fil terseancarray tsing al moving@average #f biter 
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3.12 SUBROUTINE: UBIAS 


FUR POSE: 


METHOD: 


to calculate the coefficients of a model of 
given order, choosing the best time delay over 
the allowed range using an unbiased estimator 
fiolespased on Slutsky's “theorem. The flow 
diagrams for the UBIAS module are figures 


S-2l, S-ce and 3-23. 


unbiased estimates of a least square solution 
can be obtained by subtracting the covariance 
matrix and vector, of the residuals between 
the process and model outputs, from the least 
squares matrix and vector calculations for the 


estimated parameters 
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where: 1) y is the output sequence 
2) X is the input/output matrix 
3) 9 is the vector of estimated 


parameters 


he 





gre 








» 
thers 


v 












fn I M4! 


fetbew a to esnetot iets oad sielvolao oF 


home 
“evo \eleb omtid tasd sit geetepdo, 19b10 nevegey >) 


Ji ‘ eat Ala 

todemiaedsd beretons as goatee aynan bewolia ef2 - ; 
, 4 

wort? sadT one yreesd e'ydeter® ao beesed [OF] 7 


. ' i 
14: sts sfiubeg 2fTO «4d wok semtgsis . = 


- 
Sine bars SHE ,tS=E 


Poo Ol 










: ‘> #269 
\ ; 


eee vi 
on 
4 
rl » 
7 ; wee 
of fs Wi 


4 


ou 


notiules aeavge sens! « to soteridee Bete 


~ 


‘ ‘ ae eee 
‘yavoo edt pottosywtiees ya beh ieiee Se Bee ae 

| Me 
negwded elewbiess stg ww ,Abtosy tae aithde@ | y 


) 7 7 - 
48h add mot? ,ettqauo, Esham. nate sepaotg, Ste - - 
2 es tay i. 









cto AcT anoiseswoles Tosoev bas SITTER, & 


erotome 18q be bse A 
j a ais Te 


al J) 






; _— T 

(T.#9 ; {9 oa x21 $3 - Ss 
 aoWeupes tqivo ead. at : 

bien fuqivo\syan: odd g it | 


_ - EAE: 
botentgee Te volose et mR 





143 
4H) Ro 9is the covariance matrix of the 
residuals 


5) Q is the covariance vector of the 


residuals 


J) ah = 
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Figure 3-21 
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VARIABLES: 


COM 


COV 


a real square matrix used to hold the modified 


COVM matrix 


were diueivectoreused to hold) the modified COVY 


vector 


a real vector containing the covariance vector of 


the residuals 


a’real= Square ‘matrix “containing “the” covariance 


matrix of the residuals 
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SUBROUTINES CALLED: 


BCOM P 


CCOMP 


PREDS 


PRINT 


RCURG 


RUNAV 


computes the covariance matrix of the residuals 


computes the covariance vector of the residuals 


computes the deterministic predictions of the 


discrete model 


prints out repetitive headings 


computes the parameters of an autoregressive 


model by a recursive procedure 


filters an array using a moving average filter 
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3.13 SUBROUTINE: DIAGN 


PURPOSE: 


METHOD: 


to take a model of given order and perform 
various tests to determine whether the model 
TS) ar satvisracvory fit. to the data”-and whether 
the model is a significant improvement over 
lower order models. The flow diagram for the 


DIAGN module is figure 3-24. 


the following statistical Tests are used 

to determine 3 Re the CUE ct. model ifs 

satisfactory: 

1) F=ratio determines if the mean square error 
of the current model is significantly less 
than that of lower order models 

aby beter minant of the covariance matrix, if 
zero, indicates a failure in the inversion 
of the covariance matrix (the model order 
may be too high) 

3) Volume of the covariance ellipsoid gives an 
estimate of the volume, in parameter space, 
in which the true parameters are likely to 


be away from the estimated parameters 
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Parameter significances are used to 
determine sae any parameters can be 
eliminated 
Autocorrelation function of the residuals 
determines whether the residuals are white 
noise 
Cross correlation function of the input and 
the residuals indicates whether the 
residuals are independent of the input. 
Residuals should be dependent on the 
disturbances only, there should be no 
dependence on the input if the time delay 
mancorrecat. 
Cross correlation function of the input and 
the prediction errors indicates whether the 
prediction errors are independent of the 
PID UG. Fredictton “errors should be 
dependent on the disturbance only, there 
should be no dependence on the input if the 


model order is correct. 
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VARIABLES: 


ACF 


ACFSA 


AF 


CCFUP 


CCFUR 


Dis 'T 


FRAT 


ISAI 


a real vector that contains the autocorrelation 


function of the residuals 


Seo teat Variavle that “holds “the normalization 
Peccvor es fOr. tne autocorrelation function of the 


residuals 


a real variable that contains the normalization 


MaCcbOre bor che Cross correlation functions 


a real vector that contains the cross correlation 


Pounccror or tne input and prediction errors 


ea real vector that holds the cross correlation 


function of the input and residuals 


aereat Variable that contains the 95% confidence 


Lam rt 


a real variable that contains a comparison value 
for the F-ratio test to determine if the present 


model is significantly better than the previous 


model 


an integer vector that holds the position of the 


points in the cross correlation function of the 
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input and the prediction errors that are outside 


the 95% confidence limit 


an integer vector’ thate holds» thes position iofi.the 
points in the autocorrelation function that are 


outside the 95% confidence limit 


anrinteger*vector .that«holds the: position of the 
Poticseeinesunes Cross, correlation shuncvti1on .of, une 
input and the residuals that are outside the 95% 


confidence limit 


an integer variable that counts the number of 
Police cne autocorrelation function “oucsivde 


the 95% confidence limit 


an integer variable that counts the number of 
points in the cross correlation function Of une 
Papuieeeand the residuals outside tne 4952 


confidence limit 


an integer variable that counts the number of 
pomnts in the cross correlation function of the 


input and the prediction errors outside the 95% 


confidence limit 


a real variable that contains the F-ratio value 
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a real variable that contains the volume of the 


covariance confidence ellipsoid 


reablevariables used for-the calculations of 


the standard deviation of the A and B parameters 


a real variable that contains the mean square 


error of the previous model order 
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SUBROUTINES CALLED: 
DGSW Supplies messages and program flow control at 


decision points in the program 


GRAPH allows the user to plot on a hard copy terminal, 


a graph of a function supplied by the user 
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3.14 SUBROUTINE: MDISP 


PURPOSE: 


METHOD: 


to allow the user to eliminate insignificant 
parameters from the model and to determine the 
effect these eliminations have on the impulse 
response of the model. The flow diagram for 


the MDISP module is figure 3-25, 


the module has the following abilities 

1) elimination of insignificant parameters 

2) calculation of the model's impulse response 

3) computation of a residual sequence using 
the best parameters and the modified 
input/output data. An autoregressive model 
of successive orders is fitted to this 
sequence. The best model is used to filter 
the residuals to produce an uncorrelated 
sequence which is an estimate of the 
sequence driving the system's noise model. 

4) the best parameters are used to calculate 
the predictions of the output and the 
error sequence. An autoregressive model of 


successive orders is fitted to the error 
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Sequence which is the estimated output of 


the system's noise model. 


ry 
if 
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Figure 3-25 Flow Diagram for MDISP 
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VARIABLES: 


EEM 


INT 


TROMAX 


JFRSTM 


OUTT 


SIGMM 


SIGS 


a real vector used to store intermediate 
calculations of the residuals and the prediction 


errors 


an integer variable that contains the difference 
between the number of times the input and output 
sequences were differenced in the subroutine 


DMOVE {3.7} 


an integer variable that contains the 
maximum autoregressive model order for the 


minor subroutine RCURG 


an integer variable used to store and 


change IFRST without destroying IFRST 


a real vector used to hold the model's impulse 


response 


a real variable that contains the mean square 
walties of the residuals or of the prediction 


errors returned by the minor subroutine RCURG 


a real variable that contains the mean square 


value of the prediction errors for a zeroth order 


autoregressive model 
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SUBROUTINES CALLED: 


DGSW 


FREAD 


GRAPH 


PREDS 


PRINT 


PULSE 


RCURG 


RESID 


RUNAV 


Supplies messages and program flow control at 


decision points in the program 


allows esthe,;user+\to.use-,free. format . input, .and 


provides certain error checks for the input [21] 


allows «the user to plot on a hard copy terminal, 


a graph of a function supplied by the user 


computes the deterministic predictions of the 


discrete model 


prints out repetitive headings 


calculates the impulse response of a given model 


computes the parameters of an autoregressive 


model by a recursive procedure 


computes the residuals of the estimated model 


filters an array using a moving average Silver 
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3.15 SUBROUTINE: NOISE 


PURPOSE: 


METHOD: 


to calculate the coefficients of a given noise 
model using the least Squares method. The 
flow diagram for the NOISE module is figure 


3-26. 


this program produces a noise model 


-1 § 
s Pe CZaee )s (3.8) 


—— —— 


ianerr on 





Where: 1) € is the prediction error sequence 
derived by MDISP {3.14} 
2 ) is isSem ithe (uncorrelated smani yin es 
sequence derived by MDISP {3.14} 
Various model orders from 1 to the system's 
model order are used. The mean square error 
of the noise model, the noise model parameters 
and their standard deviations and the 


autocorrelation function of the noise model's 


residuals are calculated. The user then picks 
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the best noise model order and this is then 


passed to the subroutine FINAL tL Ostse 
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Flow Diagram for NOISE 


Figure 3-26 
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VARIABLES: 
NN an integer variable that contains the noise model 


order to be determined 


SIGMP a real variable that contains the mean square 


error of the zeroth order noise model 
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SUBROUTINES CALLED: 


DGSW 


FREAD 


GRAPH 


MCOM P 


RESID 


Supplies messages and program flow control at 


decision points in the program 


allows = =tirem user’ "to ®uise’* free =format~ input = and 


pEovruececertain error*®checks*for=the input tet) 


AmLOWwsm one weer 'vtO=pDLOt*Oon~a "hardcopy terminar, 


a graph of a function supplied by the user 


computes the covariance matrix, the covariance 


vector and the determinant of the covariance 


macrix 


computes the residuals of the estimated model 
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3.16 SUBROUTINE: FINAL 


PURPOSE: 


METHOD: 


LO wOoUusSDLayestne “sy stem andeenodse transfer 
functions and to calculate the zeros and the 
poles of the system and the noise models. 
This program also allows the user to take the 
Peremotsplay ODLLON called SULOR iS wus he 
flow diagrams for the FINAL module is figure 


3-27. 


the zeros and poles of the system and the 
noise models are calculated by the subroutine 
POLTRe which is part of the IBM 360 Seientific 


Subroutine Package [19]. 
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VARIABLES: 


AA, BB, 


BLANK, 


cc, DD real vectors used to store the system 
and the noise transfer function coefficients for 


determining their poles and their zeros 


BLANKB integer vectors used in object time 
formatting of the transfer function to eliminate 


unnecessary numerical signs 


BRAC, BRAC1 integer VGC,0O1.5 6 Us © lamer) object time 


formatting to provide a bracket at the end of the 


denominator 


FMTA, FMTB, FMTC, FMTD, FMTAA, FMTCC integer vectors 


that contain the original formatting of the 


transfer functions 


LER an integer variable that is the return error code 
of the subroutine POLRT [19] 

JUN1, JUN2, JUN3, JUN4, JUN11, JUN33 integer vectors 
that are the work areas for the formatting of the 
transfer functions 

LINE an integer vector used to store the format of the 


line between the numerator and the denominator 
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MINUS an integer vector used to store the format of the 


exponents of the transfer function 


PLUS, PLUSB integer vectors used in object time 
formatting to provide a plus sign for use in the 


transfer functions 


ROOT1 a real vector used to store the imaginary roots 


of the polynomial returned by the subroutine 


POLRT [19] 


SAVE a real vector used to store the real roots of the 


polynomial returned by the subroutine POLRT [19] 


SAVPOL a real vector used to contain the absolute value 
of the polynomial coefficients to be used in 


formatting the transfer functions 


WORK a real vector used as a working area by the 


subroutine POLRT [19] 
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SUBROUTINES CALLED: 
DGSW Supplies messages and program flow control at 


decision points in the program 


BUGaeets. lb supplies, eraphical displays of data files 


POLTR computes the real and the complex roots of a real 


polynomial [19] 
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3.17 SUBROUTINE: FDISP 


PURTOSE ¢ 


DATA FILE 


1 


allows graphical display of data files and 
some file modification. Gives the user the 
Gpelon of major filer-modification, by linking 
LOsmenbe subroutine FDIS2 {3.18}. The flow 


diagram for the FDISP module is figure 3-28. 


NUMBERS: 

U, unmodified input 

Y, unmodified output 

E, errors after DAGEN {3.6} 
U1, modified input 

Y1, modified output 

E1, residuals 

a predictions of the output 


E2, prediction errors 
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Figure 3-28 Flow Diagram for FDISP 
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VARIABLES: 


AA 


ALPH 


BLANKS 


1G 30 Se 


DPF 


IF X 


INIT 


a 


a real variable that contains the mean of the 


file being considered 


areal. vectornmithat., contains. the. tittle ;osthe 
abscissa designator and the ordinate designator 


of the plot 


a real variable used to clear the real 


vector ALPH 


an integer variable that contains the number of 
differencings of the file to be done in the 


subroutine FDIS2 {3.18} 


an integer vector. used to store the file numbers 


TOD eG. DLOLCed 


an integer variable that is used as a switch to 


indicate if a particular plot is finished 


an integer variable that is used to determine the 


type of terminal that the plotting will be done 


at 


an integer variable that is used to calculate the 


time constant (s-space) for the high pass filter 
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of the file to be operated on in the subroutine 


PUI s2 {3.18} 


an integer variable that indicates the type of 


graph to be plotted by the subroutine CGPL [20] 


an integer variable that is used as a dummy 


variable in the subroutine CGPL [20] 


an integer variable that contains the current 


number of graphs on the plot under construction 


an integer variable that contains the maximum 


number of graphs to be done on the plot 


an integer variable that contains the order of 
the autoregressive model of the file operated on 


Pimecne subroutine Foloe 13.183 


an integer variable that contains the number of 


files to be plotted 


a real vector containing the autoregressive 


parameters of the file being considered 
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X9, Y9, Z9 real vectors used in the subroutine CGPL 


b20) toe storeisthe-value. of..the sabscissa, the 


ordinate and a third parameter which is a dummy 


variable in this program 


“ie 


in 
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SUBROUTINES CALLED: 


CGPL 


DGSW 


FDIS2 


GRAPH 


plotting routine supplied by MTS applications 


subroutines [20] 


Supplies messages and program flow control at 


decision points in the program 


13.183.-does*-file. modifications and calculations 


on the current file 


allows the user to plot on a hard copy terminal, 


a graph of a function supplied by the user 
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Be tomoVUBROUTINE: FDIS2 


PURPOSE: 


METHOD: 


to modify, to find auto and cross correlation 
functions and to fit an autoregressive model 
tosthe file being . considered. The flow 


diagram for the FDIS2 module is figure 3-29. 


the program uses the following to modify the 
rite: 

iyenrenepass filtering 

2) differencing 

3) moving average filtering 

The program also offers plots of the auto and 


cross correlation functions. 
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Flow Diagram for FDIS2 
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VARIABLES: 


CCF 


EE 


3 


MFRST 


NU 


NOMAX 


NUMF 


a real vector that contains the cross correlation 


function 


a real vector used to store the second file for 


the cross correlation function 


an integer vector that holds the position of the 
points in the cross correlation function that are 


outside the 95% confidence limit 


an integer variable that is used to contain the 


number of the first data element to be used in 


Picenaghapass filter 


an integer variable that contains the order of 


the high pass filter 


an integer variable that contains the maximum 
autoregressive order for the minor subroutine 


RCURG 


an integer variable that contains the number of 
the file to be used in the cross correlation 


function calculation 
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a real variable that contains the mean square 
value of the file returned by the minor 


subroutine RCURG 


a real variable that contains the time constant 


ofe thnethigh pass filter for the? file 







- 0 
: iv 
ORV OR. Jonnie ‘Ore 
..o8 2 
Oh 


mane an 
- 





Pits “at BE Sty wise saad oldstasv | 
he 4 af, 


Stee and sey *oH227 sued Agta’ wie 


5 


149 


SUBROUTINES CALLED: 


DGSW 


FREAD 


GRAPH 


RCURG 


RUNAV 


Supplies messages and program flow control at 


decision points in the program 


allows the user to use free format input and 


provides certain error checks for the input [21] 


allows the user to plot on a hard copy terminal, 


a graph of a function supplied by the user 


computes the parameters of an autoregressive 


model by a recursive procedure 


filters an array uSing a moving average filter 
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CHAPTER FOUR 


PARAMETER ESTIMATION OF A SIMULATED PROCESS USING UNALIP 


4.1 INTRODUCTION 


11s section compares the performance and 
BaemecOmputratLiTonal timle of .thé four identification 
methods used in the UNiversity of ALberta Identification 
Package (UNALIP). All the methods are off-line 
identification routines that lead to parametric models. 
This chapter has been organized in a manner similar to 
the formats used by Isermann [2] and Saridis [28]. 
Simulated processes were used because: 
1) the systems and disturbances were known 
exactly and could be duplicated 
Pomtne. inpute and" the, notse/signal ratio ctype 
could be varied easily 
The methods of identification used were: 
1) Generalized Least Squares (STIME) {3.6} 
2) Maximum Likelihood (ASTRM) {3.7} 
3) Instrumental Variables (IVOFF) {3.8} 


4) Unbiased Estimator (UBIAS) {3.9} 
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4.2 TEST CASES 


Three simulated, linear, time invariant processes 
were used as test cases. These processes were: 


1) Second order oscillating process 


1 


2, CoN, iia 
ria) 
1 x 
= LeU Zz, +O ee (4.1) 
ya pete Fae ay ame 


2) Second order nonminimum phase process 
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3) Third order low pass delay process 


ath en 
G,(27') fey a Bg Sr") 
fiteaes) 
= = afi 
A 0.065z + 0.0482 - 0.008z (4.3) 
a ee Wl = os kor 


Figure 4-1 is a general block diagram of the simulated 


processes. 
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Figure 4-1 General Block Diagram of Simulated Processes 
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The output, v, of the system to be identified was 
contaminated by an autocorrelated discrete noise 


sequence which was generated by a noise filter Grice 


).. 
The filter was driven by a discrete white noise 
sequence, w, with normal distribution, zero mean and 
variable variance, “ftp depending on the type of the 


noise/signal ratio desired. The same discrete transfer 


function used by Eykhoff [9] was employed for all three 


systems: 
= 1 + mien Ss 
GeGzente. “She rete 
* -1 
os Leal, 
Z 1 (4.4) 
1a ALIe te. a Ofenas 


Maes discrete white noise input of the noise filter 
eam: was produced by an IBM System Subroutine Package 
[19] routine called GAUSS which generated Gaussian 
distributed numbers with zero mean and variance as 
specified. Therefore the noise sequence was 


reproducible. 
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4.3 IDENTIFICATION 


The input signal to all processes, a pseudo random 
binary sequence, was generated in DAGEN {3.6} ee 
accessing the minor subroutine PRBS. The pseudo random 
Dananry sequence, thegraweinput>—had=-a=ength’ of 511 and 
the two states produced were -1.0 and 1.0. Figure 4-2 
is a plot of the raw input sequence used in all runs. 
DAGEN {3.6} took the raw input sequence and passed it 
through a specified process model which produced a raw 
output sequence. The noise sequence was generated and 
added to the raw output sequence. Filtering of the raw 
BADUL OY Output wassnotsutilized anethis study. The raw 
input and output values were passed to DMOVE {3.7} where 
Pne values were transformed. into._.-mean. free modified 
input and output sequences. No. tritvering -ot apne 
modified input or output was employed so the modified 
input and output values were used by the estimation 
routines to determine the process model. 

Tests were performed for two noise/signal ratio 
types at input/output sequence lengths of 100, 300, ra 
500 data pairs. The error of the estimated parameters 
and the error of the resulting impulse response for the 
estimated systems were tabulated for all runs. The 
errors were calculated from formulas provided by 


Isermann [2]: 
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Figure 4-2 
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128.0 192.0 256.0 320.0 4.0 448.0 


64.0 


SAMPLING INTERVAL 


Pseudo Random Binary Sequence Used as Input 
for Simulation Runs 
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a) Squared relative parameter errors 


2 11/2 
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where: 1) 554 is the squared relative 
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b) Mean squared parameter errors related to the 


mean squared true parameters 
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where és is the mean square parameter error 
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related to the mean squared true 


parameters 
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c) Mean squared impulse response errors related 


to the mean squared true impulse function 


values 
eae 1/2 T ee 
Agy re age 
t=0 
SS = ee = IE ea ORS Cae 
Ts 2 
s 
Bt zr By 
t=0 
where: 1) Ag, = B. - & 


24) Sy is the impulse response sequence 
of the true system 

3) g. is the impulse response sequence 
of the estimated system 


4) T LS Lnes sett Ling. |times so temtne 


s 
system 
5) S5 is the mean squared impulse 


response error related to the mean 

squared true impulse function values 

To obtain a better comparison of Ghe four 
identification methods, five runs with different noise 
sequences were used for each noise/signal ratio type. 
The noise/signal ratio types used were low and high. 
Low or high noise/signal ratios types are generated by 
manipulating the standard deviations of the Gaussian 


Mesgeraobucion driving the nolse filter: Noise/signal 
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ratios are calculated using the formulation originally 


encoded in OXIP: 


Vid] (4.8) 


The noise sequences were generated using a different 
seed for the subroutine GAUSS [19]. The seeds used were 
1, 3; 5, 7? and 9. The noise/signal-ratio values for the 
five different noise sequences and the three different 
processes are given in Tables 4.1 and 4.2. The standard 
Ggeviation evalues for the .low.and high noise/signal 
BaclLosw were e0.t1 and 5.08#for wproaess. I and 0.01 and 0.1 
for processes II and III. The reduction of the standard 
deviation values for processes II and III was due to the 
high noise/signal ratios generated using the standard 
Pevinution vaiues O.1 and 5°70 for®these processes. The 
standard deviation values were reduced until the 
'noise/signal ratios were similar fortes three 
processes. 

To obtain results that could be compared to 
Isermann's [2], the standard deviations of the errors 
were calculated using the following relationships for 


the five noise sequences 
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The average standard deviations of the parameters were 
also calculated and used as error measures uSing the 


expression 
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Each of the parameter identification methods, STIME 
3-58, MASTRM 43.7}, FVOFFi-(308} and UBIAS 43.9} have 
been applied to the three processes using: 

1) two noise/signal ratios 

2) five different noise sequences 

3) three different input/output sequence lengths 
Therefore, 4 x 3 x 2 x 5 x 3 = 360 identification runs 
were used to compile the results. Tables 4.3 to 4.6 are 
tabulations of the standard deviations of the error 
measures for the four estimation routines. Tables 4.7 
to 4.10 dare’ tabulations of the average standard 
deviations of the process parameters for the four 
estimation routines. 

The .computational ‘time*® of the identification 
routines was roughly estimated by determining the amount 
of CPU. time (seconds) required for an identification 
method to estimate the parameters of a particular system 
with a fixed number of data pairs. THis est imaricone.s 
similar to the comparison done by Saridis [28]. Two 
cases were used to obtain data for computational expense 
for the four estimation routines 

1) low noise/signal ratio and an input of 100 
data pairs 
2) high noise/signal ratio and an input of 500 
data pairs 
The input data was the same as the other identification 


runs. The seed of the noise generator was 1 for all 


ad 


—~ 


~ 


19 Jug@ag ap Sas ‘phys ie 





Juqnt-. nB bas, Otay taain aii 


7 


oosahanntioas 42 Paty co. 


‘anfttéoe@bht Gab = Se @ x & € x ,sq0% 


‘tcigiaq # Yo 91sdemsrx9q odd os nbewel'el 
f pees 


S$) 2ibissé ya sacd -nodtibaase #8 


snoijdatugeon woo £268 aisrdo ot 






















| | aa z 
og sordd of9 OF Bett 
os =i st " h 
pt is agte\aekoa owe 


m bie 
i 


s 
a 


oTaw me fone sueseR teh 


reupes ug s4oNSullie: Sue18htth son 


i ' ie 
esianT .ecigees eft efiqmes ot OBey 
; 3 7 _ 4 : a) i _* j 
wRhod’) Hi vah bile # irs one (Yo, anon gate 


iy 


£ oy: %+ ~ dao? fart sea wot ‘ae bial ae 11 


esigva oa4 Fo shot tetased ana 


bo} ws ans 789 BSeoo% ¢ aude Page Est 
iis i’ a fee? 


err 56 my rn 

a¢ 
Yoelp if Sao Io -eai? choke atugaos Be, id 
J=5 4S re ssalsee aouatle ics gis tis 


ot ne. aed ber tepes (abnovss vs et nd 


y 
ib ‘eh 
6 


ce 


ze «slay 32280 shaded ie) — ieee 
ip Sigi 


uu 


; pravtias soStam 


=} 


on 
- 






163 


6€S°0 12S ~ 6LS°0 


| 
| 
| 
| 98E°t 
| 
| 


| | | | 

| | | | 
ofesg | GlGss | PLE*s LG6G°h | 9z6°2 ! 0 

| | 

| | | | 


OLX66L1°L OLXLLL°L OLXLhth Lt OLXE66°L 


WE Oy lomo ae eee eine BEER Eine PE eesemelinesl AE at scerenay| cr orf 2 oon 


og 


OLXLOL’S OLXQS9°6 


| | | 

_OLXz9on'L | OLL*O 6LL°0 OLXLEL‘L _OLX6g0°E -ObeLHR?S § - “of 

z | abe € | Lz | 
| | | 

| | | 

| | 


LOL‘O 0n2°0 992°0 OLxXgLn’s _OLXL9S°6 _OLXG69°1 ¥ 


Gat 


3 


® 


8SE°0 Lott 


| 
| 
| 
| 
| 
| 
[= 
| 
| 
| 
| 
| 
| OLXLQn’S 
| ------------ | ------------ | ------------ 
| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 


| 
| 
| 
| 
| 
E | 
A ata yon Vp 
| 
| 
| 
| 
| 


| 

| 

| 

| 

| 

| 

| 

| | 
| | 
| | 
| II ssoeooud | 
| | 
OLXOGE'L 2 OLxzge'2 | 
| | 
| | 
| | 
| | 
| | 
| 


5 Ran pean SEP ag S end TE Ph dchebitetksths te aptly Geter re 
| | | 2%, | 
6E£°O | Cth’ O | 699°0 _-OLXOL9'E _-OL¥She°9 va _OLX6LE‘L | 0 | 
L 3 
g60°tr | 68r°l | 86L°2 ,-OLXBLE TL -ObL¥LGZ*2 - _OLXZ9L"*t 7 % ) I sseaooud 
| | 3 
San) | LLh'O | 69S°0 | _-0b*200° € | ee _OLXZ299't | prob XO8E*L | Zo | | 
00S | OO€ | OOL | 00S | OO€ | OOL | yySueT und | 
yuSty | MOT | TeuSts/sestou | 
~----- ee ne a nn nn ee ee === | 
AWILS | 


(AWILS) Souenbs 4ysea] paeztTeueuey 
JOJ seunseoy uouuyq ey44 Jo suotyetasqg puepueys ey4 jo uostuedwog €*°h ATAVL 


















lt alactnieaeiachaeaieal cares 


ee 


‘ | 
} “Oreo! .2 
i 


~ a 
q i 
—. oa 
r eae .e 
—— 
~ 
aa @2& = ae @aownw ee Oo ee 
“ve 
‘. 
so «@ 
Lar 


a 
peice wate moe [- aoe nee ne 
€4e. § 4 


ter .¢ 
eae.¢ 


TPBS0 


C2 +2 


> 


2. | 


‘as 


in die nik apenieesy 


out 


nawoaw Wao e @ 


oe 


>! 


- oa a 


us 


oo 


! 
' 
: 
) 
; 


e 


4 


teqee~ x 





“OrxS00.£ 


“Orxr ee 


“orxt ye 


os #BEE.1 
Or earRve 


* 
—-— eo = Oe Or eo = 


q 


* * 


anc tgp et lila ee pe ee 





¥ 
antec} +ys == 


; 





be. 
~~ 
@ 


_— 


“Gr 81S, 


ud 


ee - seeowe 


“Stacee- 


4 


“orstée 


OF neao.€ 


yteted, 


> 
——* 
“ 
~ 
: 3 





or xt rer 7 ; hae 
Laas 


rT erasae.S | 


ai 


“Orxbea.' { . gf 


a un atime 
a 
— 
- 


- > acme es =e - stigens: ¢ >< ae Pe 
raTee.t erxee?e.t |. ,0 | 
| & 
22.0 | asevs: fi)" 
: Y ay 
: é 4 Se€.9 | 5° ; 
f i ae ioe . a 
a a a += aston kann tnee ses Red ae 
i) a 
’ ‘! - 





. 





, 


a 
sae 


164 


| | | | | [er 


61h°O tite * 0 Lhe °O 


Let 


5 


| 
| | | | | | 
| | | | | | 
LLo°s | yet “SG 7 fH96°L | Chobe i G6L°E€ | © | III sseooud 
| | | | | | 
_OLXH L6H o1x000's | Dot Of ce OL XO099 6: 41 Se OLXE nT: 4 OLxEe2z°E€ | | 


Zz 


20326 eS ae |------------ |-4---------- |-4---------- [-2oSeeee SSeS |------ |------------- 


| on 
O1KEQE* _OLXLHZ*°2 | _OLXZL9°8 _OLX60h*L ¥ 6 
z Zz 
| 
| 
| 


| | 
| | _OLXSEn’? _OLXLE2?*H 
oe lees 
| | 
| | 
| | 


: Lg 


| 

: 
_OLX9EL*S _OLX6SL°9 LOL‘O _OLX6S1L°S Laue? 
| 
| 


z 4 4 


3 


¥ 5 


| 
| 
OLX IG TO LxGh OA 
| 
| 


OLXglS°9 90L°O 88lL°oO 0Lx980°9 OLXSOO'L OLXLS6°L 


| 
| 
| 
| 
| 
| 
s 
| 
| 
| 
| 
| 
| 
fe Se ee a aaa 
| 
| 
| 
| 
| 
i 
| 
| 
| 
| 
| 


| 

| 

| 

| 

| 

| 

| 

| | 
| | 
| | 
o 7 TI sseooud | 
| | 
| | 
| 

| 

| 

| 

| 

| 
| 


On2'0 ! 6SE°0 | 955°0 | _-ObxHSHE | _-0L¥1L88 °S | OLX LONE a + | 
618°0 | ger‘ | th6°L | ,-OLXELZ I | ,-OLXES0°2 | ,-OLXELL'S ! mi | I sseooud 
eee eho} Is ebro 1s cleo’ |, 0S | ee 
00S | OO€ | OOL | 00S | OO€ | OOL | yyZueT una | 
Mena stg oe A ee ee | 
i eee WULSY : 


(WHULSV) POOUTTOEATT unutxey 
JOJ sounseepW Jouug oyy JO suoTyeTAog puepueys oy Jo uosTtuedwo) i'h ATAVL 




















cart aps mop snenessee eee tee 

















_ MAEREE 





OFx€is.t 





ee eg PEPE ® cuca — wer ee & eh =orgerer we ~nbeatiat =n saeptePA See ° 


ithe 





Vue Gahws ea ee a, — ee — oe 


mires ar 4 


if? 


ee 


——_ oe a a ee ee ee ee me ee 


ak 


r 





t “Oted ee. i ] 
i 2 
“Create. 4 % 
+ i $s 
“opetes.# 1. 
| st q ae 4 


-— =: ee ai mes 


er ety I re teen: pee 


“ohetSSi€ dye | lige 
Lao pm 
“OT .€ . ie » ss 
| a an Se 
ree.a°q a 
P — a a ; 
ewe ee ee ee ee ee a ee eee oe 


*~ 


165 


| | aes Pd 
9€t° 0 60S°0 | re) 
| LZ, 
Fisee ss 
| 


| 
| 
Lege | LEE 
| 
! 9 


| 
| 
o | III sseooud 
| 
| 


© 
foe) 
iat) 
— 


_OLX960°S o 


pHa eg gn ae pe 


_OLXt61°8 ~OLX9tO'L “OES Gcc cc 


| 
| 
| 
| 
| 
~-f--- eon nn | === Se ee an em ew mi i sm 
| 
| 
| 
| 
| 
| 


CAR 


9 


POPAL ES Tf SOPXTSE c _OLXS86°t ~ 


€ € 


,- OE *SS0°S es0l cau g ,-Obxech | 9 


3 


* 6 


| | 

eee 
ree 

hj at et 

| | 

| | 


ge Ol*GRE 9 phe Po ae: ge O Lie ToS 


| | 
| | 
| | 
| | 
| | 
| | 
en ce a a ge 
| | 
| | 
| | 
| | 
| | 


| 
| 
| 
| 
| 
| 
| 
| 
| 
II ssoooud | 
| 
| 
| 
| 
| 
| 
| 
| 
| 


ELL*O g62°0 | g8S°0 OLX Onn € -OLXZE6°S | _OLxgge't | ar | 
LE9°0 r60°L | bht*2 | | _oLxosert ,-OlxGhh*2 | OLxgls) | a | I ssoooud 
Bp eB £9820. |g SONsocere | yCLiDSO.d. | plete t | es 
00S | OO€ | OOL | 00S | OO€ | OOL | yyZueT una | 
“REE iia aaa emeiiaeaieniiianay Siem tra, ns. | 


(44OAL) SOTQETUeA TeJUSeuNAySUT 
JOJ sounseop Jouuq 044 Jo suotyetasg puepueqys ayy jo uostuedwmoDg :G*h ATAVL 








iSO | 








Peay 
62s.9 


“sore 


FED 


aa 


r= , 
“OPxSee .s 





be J 
“ee & 
ie Bl. 


ere ae oS eee ——— 
enmrew i 
: i “Ores C..< ! Or eee 9 
f : ‘ ’ 
see.k | vet. 6 
} ; } 
ost.@ | ebe.0 | 
{ 





ee art rarer a ar 


eer.0 


-. se ee ee ee) 
4 


“GP xosd,, 2 


—— a. a 

tes 3 7 
i t. 
; U << e 
" ° 
t Wax)? 
e! 
: 
‘ 

i « . 

{ . 


Sep ~—-- ——— 


Ta ReSE- ¢ 


eee “8 =e ee 


ae A 
Gixert.! 


“@h2h28. 0: 


ree. S 


oo —_ 


ch ee os Ml a te ay ts Ol A CE 





‘~gpeeto. it 


t 


y= wee eo eS & —" oe eer 


“gtaess.S 1 


— ——, <<, 


166 


svian | 


| | | | | | | 2%, | | 
| L29°0 : 109°0 | €69°0 | LOoS‘o | 6LG°O | H8SG°O | Kinet | 
| | poh | | 
! 2SH°S | nhe’s | 668° t | 98E°h | VeGet | L96°t | ar | III sseooud | 
g 

| OLxz60°S | ,-OLXLSL "8 | OnlL°Oo | - OBXL Lei | ObX kre be} ObX ELS tam | ae | | 
| -2---------- | -2---------- |------------ |-4---------- |-42---------|- 4---------- |------ |-------------| 
| | | | | | [jez oo | | 
PPelOLXEGEStO IT “olxechscal? _OLxh6gnsel) _o1xcséstel! _olxnhsscal} olxgooesa! o | | 
Ea ip [?* [Bie ee bes | Lie) | 
[P= 01x00t- Go|] =-01%669" 9-4 066° 0-41 £-=01%250+G6=-| a-01%0092e-- fe 04x9z)=1—4 9 | II ssaooud | 
a (a | Be Ris i | a | | 
FP _OLXH88°9 | ELieo | th2° 0 be -OUXHLLe9 ‘= _OLX0Z0°'L he _OLX96E°Z | 0 ! | 
| | | | [ | Be roe J | 
| ELbeo~ coero~{ Les‘o | _-ObXnEnE | OL XZL6°S | OFX 982" 1 | o | : 
| | | | | le ke! | 

! 2S59°0 | GLL‘L | cples : i _OLXLH2°L | = _OLX9EL*2 | fe _OLXLEQ‘H | ret: ) I sseooud 7 

g 

| 89L°0 : L9z2°0 | L2n'0 fe -OLxUUeeee ! -ObeG0g 5S | a _OLXZnE*L | 0 | | 
fe 00S | OO€ | OOL | 00S | OO€ | OOL | yysuetT ung | 
| ustTyY | MOT | TeusStssestou | 
| 
| 


(SVIGN) woyeUTISY poesetquy 
JOJ seunseoWq JOuuq 944 JO suoTyeTASqd puepueysS syq JO uostuedwog :9°t AIAVL 








a : —< “oe 


“POS.0-f. 7See 
Se ee ee ge 3 
: eH: i: ett is || 24 

: “4 ‘ : y & " C3 t 
; 

} 

i 





tee 





€0€.,90 SF re? .0- 
a ; “a * 3 - 


eqroewteernw es & fom me me we se ae Sainetin te ‘caheieied eo 





. 
2-2 «&— 4&4 & & & 


= - 
j : <yoen-- = = 


£tt2-] (2220 | Of«<eT?.a Srzor9 ./ 


t 
“or mbites + Pe 
+ BABS. Peaitel | 9 % 
| 
. 
| 
| 


- 
bd "y ‘ «a 








; 4 ae 2 —_ 
2 Tenner, * 4 oy .6. Ff “ovyetet.2 | “orxose.4: 
. 
, 


; “Oreesi.s Of .a reece Teh TE.S 5 er . —_ 
: eS a sy ee 
SS OEE nO Ae ae ee ets gers sar iers> | 
t  GhaSeO.e° | STeVete i SPP o | DFst* | eisver. tT Orartre 47 4 _—s i, 
2 i : . | P = : 
‘ J 


' 
6 
: 
” 
wy 
ae 
we 
w 
i 
<> 
“a 
~ 
& 
‘© 
= 
ee ee ee 
" 
- 
“e 
— os 
el en, te. 
a 
i 
y 


ep - 6 o ESS.O G 6 ere. bol 29 
aq ; | sz |. 
n> & eo & ~— +. hee = = a oa » — = - : tty eh es a a a Sa ee a eee ee ee 
aes 
amy 
“ 


167 


AWILS | 


| ,-ObX99L°€ | ,-OLXE90" tr | ,-0bxz08"€ | ,-01X682"0 _OExeteEO Ee [1 SOeRoh Sa: | | 
| 7 -OLX66L°E | A OLXLZO"N | POLXSLO"n | F -OLXGOL"E l) * OuUxGOE"S |) = ouxakicves ©] o | | 
| *_ouxigg"€ | *_orx09L°9 | f_orxono'l | f_otxnne-€ | (otxerge'9 | “_OLXELO'L | a | | 
| Loz‘o | ggLo | HlLZ°O | oo roa | hema |e 2i:"0; | 32% If EER SaeRoade | 
| EL9°O | L99°0 | 659°0 | on9°O | €99°0 | eters Wie ae | 
| €SGS°0 | €9S°0 | €2S°0 | G8n°o | Z20S°0 ! EzE'o | o | | 
maar mayen me ems | aera mg ma mc a ren | a eee | acer eee 
| orxzgi-g | |orxesis | orxzenst | |ovxmmte | | orxuge's | _orx6na-t | “eas | | 
f fsoumeess |) “_o1xese~9 | nlh°o | “_orx6z6°e€ | “_orxz9g°9 | *_oLxoor'l | se” | | 
| *_orxtta'9 | Shirin 6 | -o1xnl6°6 | *_orx6zc*t | *_orxstL-2 | “OQ ang0 GI) fo I Uk SSannae 
ii A OLXLGL "6 | HEL *O | OSL°O | pO Ex26b 2% | er O EX 69L "6 | prot ¥L89" § i] o | | 
ZI ---_----|------------ |------------ |-4---------- |-4£---------- |-4---------- |-5--- |-------------| 
| gsr o | z2g9°o0 | G66"O) |i _ounxeoevo; li -onsteois lt _onshoave | “Bing | | 
| Keo | gLe°o | o€9°0 | S_orx€gore | *_oLx96n°S | *_o1xzL6°6 | $30 | | 
| Gez'o | QgLE°O | fete 0) |, * Sane | S_ouxeas*t | S_orxe6s~s | re | I sseooud | 
| 90€°0 | Ere *0 | LG 0 | - OLXZnO0° us | Ke OLxzZgl-L | ‘a OLxLen'9 | o | | 
| 00S | OO€ | OOL | 00S | OO€ | OOL | yyZuetT una | 
wan nn nnn nn nn nn nn nn nn ne | nnn nnn nnn nnn nnn nnn errr nnn nr ncses [ner rcsccrecsssscss | 
| yusty | MOT | [TeuStssestou | 
| 


(AWILS) Souenbs yseey paeztTTeueuey 
kq peqyewtysq Suaeyemeueg ouy JO suotyetaAsqg puepueys e3euaeay oy4 Jo uostueduog :)°h ATAVL 








I 
ttmwenn [nove sen zamnd 


£8 j Tas, ) 
@ft.o |. BELLO = 
~ “884.0 °| - G22.0 
| $80.0 | £20.60 


prerie aedl Ketetetaietehtdechaetalt tasiematasehcicie tai’ 


.. GbE ort - eat a. 


Ma “cr xsie. f: 4 vateagt 
““ptxbee. ? 7: ce 
+ ORES. ate 
“OPebde | he 
operat eer eps aets 


, O1x5er, 1 “pixeaT.2 4 . Ofx 8a. 2 = 














| ,_oreter.e 


= 





. 
} 
Ot ateBee | “Rixlte.e | .“erxese.' 4 ,lOieerT.S }., OiebsO-e” 
“5 Abar bbage 4 ,. ORS |, Otedse.e Pf etses.o |] _OFRAOT 1 
Gegeee. "Stayer t 4: Otever.s j oTeTtee.c { orxges.t 
Atha bel ee Pn ae tlle « 7 | —~ Hae ee (cere iaRe, 
oe Far.o | cep | 25.0} So2.8 | “Ese 50 | 
A iige of ae | a2o0 | Gf0.0 | E400 | ve “SER bs 
. ‘oe gtr .0. | bs 0 ors.e |, s.0 |, “‘Satee f 
bs & | , Ofxeoy.a | GiadeO.t | ."Oretdelg |... ORB ee | | Of eEMRe 4 
Sew tsa snad bs | . ofytso. |, @fr87c;8 1 4 cet x , OFxeGEcE | , OFastr ss 4 
‘¢ “SPeaet.& | OraeeG ww | Mires s: | O' xO BS OPT TOL | Otatoe.t 4 
“y een 
ee sees ~> ~~ ee ee ee Se ee oe ee ee ee ee oe oe or ee ee ee ee We eae 


rw 


vai 


168 


soexsogyve (| fomonie ze || 
‘ 4 _01XS06°2 OLXLLS* 
"Torxn6z"€ | *_orx6Loth | *_OLxgLe'h | enna | -01xzG6"t | _olxLen 2 | Ef] 
_-01X006 E al "* coos ig || Z~OLXLHE'L 27 PX6E0 wt 1] “CoOPRoU Ye +| SSOUxXOGe te a] eg, | 
Bizeogis  ioeto |” te ge'to | *Torxeee"J | "Torxso"L | *louxsze-t | cee | | 
LL9°O | L2L° C - ; Oot to ahs : Ee 
asi 0 | eek | age t\* Bit 0 tl ofh “0 | 0c 20 tl eee ee ee 
------------ |------------ 08S°0 | 09L°O | 62£°O | gL5°0 | Le. | 
,-OLXSHEE | | _01x6EQ"S aera reeresteal aetertoore aeerereeati ral | 
Berta | flouxzd-9 | Florxer st | CroMictg SE | ,-01x6L6°S | olxezht | [qe | eae | 
FloLxeze*L | f_orx6or'2 | f_orxoLe: 8° ‘TorxenL’9 | f_orxtgort Ld, | 
orx9gs't | *_oLx j = bem | “sorenog-h tl Olegpes | * 2e | | 
= Sa arin one il “sovxoie ty W “xoxo || “xopx "2 | oreneceh | fe? | Ir sseooud | 
Coie *er) ee ng ¢ : oe ees | . 9 | | 
ME 00 17 GL6°0 | Si okalermtagete | Yowiywtehctc haem mie |i 3) Sein |p ten 
_OLX2GL° ; Soe 
ee aur | Se | Maa REE 4 STORE | orxnos*2 | [4° | | 
ggL°o | 92z'0 | aie io: Alam ose -OLXE0n "9 | *~orxozg°g | bio | 
eGit to” || 8gL°o | autor Tl ee ee | path Ore Ft | ©_o1xezo’9 | <p A I ssoooud | 
caaeaattcte|ecnnnenetata faeecnceetatn naten eEe S| paOlXEHETL | f_OLx69L"9 | a || 5 | 
2 ee | 001 | 00 | is ain to ee 
tetcnrnnennnennnnas [nnmnnenennnecncnnnneren sonora | uqyZuet una | 
be | ee itera! \Eatribiichahene'= Iai 
| TeuStssestou | 


| wuisy | 


kq peyeuwtysy suoye (WHULSV) POOUTTS 
2 qoweueg 944 jo suoTtTyeTtAsq puepueys te sity, iS VOSS RIES 
; 9 ?8°t FTV 





Lhaedat er ere ees 
gin 1 <a Boge tine a, a 


as 












-—— wa dknaapon tae ‘ 
Py 2 leatak “oe a ae je 
oe ‘por ce J 
: | | Hesston wena a, aR EE o\{- — 4eaz 


See | “err. s. To tere my 





one BRST. EG. Gera A ‘% <a 1 es 
Eve. o- aD 4 “praset.: ‘4 erat ott 


14> ae = Se 


. 






eat Seooecow ee — ‘> , soared ie eee: eee o 


[ 
i hi |, Orsone se 4 WE aca tage 4, Otatsa.F ¥, *98T.5 


as 
; 
xESE~f , OF eeet.s ee Leet eI - ~Oretot- ji fron sdts (E : Bite o 4 
“fpixSTT..2 sores ey « iy | (he > - , - eg ay Ps 4 2. denn ‘i 3 
= te ~ wa a hd 
ire Preset - | “Gteeer-E | ~-STARFE LS 4 ‘aes so. 
n= = & & ey > - ~—aeeaen ew = © -- |--eenseees*s ce 


saan : ee ee ee ned 


+ oe , | ‘Gie 6 =| Gor.e2 | gsc .9 
‘S *.¢ ni Sar -U j ; stot . a #. 0 


43 
4 
+ |- 
de ee 
a 


a. 


"e 
5 Sieg nine - =e 


io re 

; ee 

t . > 6 ‘ a5 
i 0 cf cxPee.y 7 . mes.0 } . orn q 

; i. i 20 i. . . *. & | fo . .- a ye ~ ce? 
t “ohxeean a | Otarer ot 7 at 46 | . Ofsero.s | of SESE Wf Bee 
To. Sreece.© | . Graevtad 4 4 Ds . oreatTPr.& 1 «Sf oO! 3688. S at Dy) 
+ “ptxiis.£ | ~o T orenee.4 | Otarhes ae” 
St Ri 


Seeger eee neem ~~ nmavwebawewnw ee . -_ a son @=a —_ P ne = _ 
; - = 


169 


AAOAL | 


| ,-OLXLLL “2 | ,-OLX8E9"2 | ,-Obxea2°€ | ,-0bXS08 "2 [it onstpegee: || — onbtasui"e | eae | | 
| ,-OFX8H8 "2 | A OLXLHSE | POLE LHL A [= pumquin's [i “omens se |) * oposite | ice | | 
| r-OLX8L6°E | Ol xeon 'L li = ouptrgo. Wl lt * orpenerc: |)“ onpeca’o) |i * obeeuon 1h | cee | | 
| L6L°0 | Weeee'oo |F * eitee’or | eeiae "op lt age lt gzz"'o? |} 5 It Tat seecoads | 
| gES°0 | 619°0 | zoL°0 | 029°0 } 955°0 | colo. Hh ee I | 
86£°0 6St° 0 tLS°0 OLn*O Zeno 68r°O | o | | 
[i i et | fm i | emcee | ren mm | imines | Epc | Fiery 
| _-OLXLE2°E | A009 oss [| oreoone tt |} _oegzeec lf _oubcogor'a; || Ober ere | pe | | 
| 01668 °E lt “Sonpcoui's IP ovpircort [1 “ ovens Ss || “oir | Sone corte | mF | | 
[Pee |S olpmeee || S_oteeter’g) |! “ oveoigtins |! “oltre |) “_oloceatios ce | 9 | II sseooud | 
} 7-0 t*995" 1 7 * 01869 "2 | 20 bx019'S | £2 0bX2SS" 1 li “Sixtus |b i oikecutye’a: | me | | 
~Z2--------- | -22--------- | -25--------- | -i---------- | $f n= | ee ono ne | enn |---| 
| goe°o | 62S°0 | ociort |] oiseentor’ge |b zorbeiehee 1) loeuaiae ot | ae | | 
| Leas 0> | LS2°0 | sogem |i So1xgenr's |} * o1lbetaie Ge [f° LoleG eer oe | ee | | 
| OL XOL2 "tt | ,-Obxn 6h” 9 | oulit’or [i “oncouter’gs Ib Sotocceee lh | f° citer Se | ee | I ssaooad | 
: FOLXEEL*E | OLX OHL 9 | Z8L°0 | "_OLXS82°g | ObXti nL | -0LXGS2°9 } o | 3 
| 00S | OO€ | OOL | 00S | OO€ | OOL | yysueT und | 
| yusTY | MOT | Teustsysestou | 
| 
| 


(AHOALI) SOTQeTAUeA Te PUSMSNAYSUT 
kq peqeutysq Ssueqoueueg 944 JO suOoTyeTAeg PuepueyS eSeusay 944 Jo uostuedwo) :6*h ATAVL 









ee 





eee st 
al dank ae’ ii epee Ne <a “wy 


; 

, Srxsee tk , OTK eet: &.+ we, o¢ ee? 
> 
; 





i "con . 
wens aaneh aoa, = ae =r i 
2 \ = oe ii Reena. Og: J ie 
‘OrxaPt. 2 str. 9..7 + OFX SBE. & 
oi. * ‘ 
2 = | 3 tar Bee ie ty om 
_BtTxade.t + ofebe2.'§ fos Stet 0. ai 
“Ot xtes. t bh te oreger.$ ne “UPeare.o 


eRe bce! aia es 
i; naewwn anak |agsnnnnn nant: ca 
» Uneeee.® |, gre 8 ¥ Woeteare.8 [az 
ane | a nck el 
~O4eeOR a 





"Ot neei.? : ee "6 Fee ht 9 Tot KETO Se 


se.€ I = Eta 58 28 
Oa 


s, 
g.,Uls 


Or xéSr. Fs : OPS O8e. 2 


TeraOes,o. |. al PT 


Si tomes aE ee wow 


<a I? o. P| 4 5 


‘ 


~~ mew ew et dt et ee eareees > ete 


é 
q~ a . os 
' ~~. V 





Ra AL I A ene ot ay le 


; 
eFa.0. |: ser.o 650.0 | 2b.9 ie Fee. 
Ts ES Og EO Pl ETS | ED gee ee 
4. OF ERD 1,  OReeO.! 2. GPwRBBik }.-q ofxfes.¢ | , OO aero -aes 
= 7 s 7 ' ’ : - ; a * = . ai t+ . 
7 : ee Wes arutec bj. . Gtater.2 4 4, 8h: .c | , Oweet.S | , Ofetesis | nage r 
oad 5 La - e "= +e ~ = } a 6 ae ~ a > . b : 
*afel yt. & l “giywhFea.s ) C F Sat : r yn < : OteTee. Ss j OTael: a | pap? } 
* ' i J - —_ o 
Se ceggsigenan claGann sme tonnn rie: a a i — a ee ee ie _———— ce em eu tn ee et ee ee ee ee ee ae ae ame 
4 Ps 
- ‘ 
: 7 
as- 
Th 2 
‘ 
7 


170 


SvIan | 


| ,-OLXLSS"€ Seco Lx cre lee olxcor st [etx co Sec ae OLX) LOeE sl er OL Oy eee eo 4 | 
| HOlxthL tt | -OLxSze*h | 9 _OLx06h hn | ~-OLXS9L"€ [MteolxsoS:e | mOLxGlo sommes cam | 
| F_OLxoL6€ | f_orxgl9°L | PAOLXOEL*L | (_orxnne € | ST OUXOL GS-9. | MeO. Oe me ce? | | 
| ogz'o | Gez'o | EQte 0 lga: 6 ics 0a ae U2 Zee ee n92°O | 539 | III sseooud | 
| z61°0 | GLL*O | 006°0 | on9°o | €99°0 | Pokxelbcses | | 
| 865°0 | 61S°0 | L19°0 | Gn’ | 20S°0 | GGG=0 mle | 
wenen nen —~—— | ---S- ~~ = =~ = | nn nm | wm mn nn = | nnn [mn nnn [enn nana | 
| ,-ObXG02°€ | _-Ob¥enS: S | _-OLxXE6E 1 | Peace tae | OLX G97Gish seOlxece tele eco | 
| _-0bXn88"€ | --OLXLGL*9 | *0LX0S0"L [Stearn 26: tl OLxe90s Ou ee X66 Ome iy: | | 
| ,-ObXT62 "1 | OL X88z 2 | = 0LX6SS°9 [eer Ox LG las tae (hee Ok Gul Late ae eine LX Ot Omen So | II sseooud | 
| F~OLXH09 "| | 7TOLXgSL*2 | FT01x198° G | £101 x65" 1 | STOLXSHL* F 2° OlxnEen’S | LD | | 
aos eneee ae hea SS6ee eee 66 [oVMasaseSaeseore ees e 6S eeeccnce eos sees Se ey SS | SSS eal eee este aS eee - | 
CLE: : ‘ : zq 
| Byte: | scone Pee at RSS he camcsere tt Paar TS | | 
" - * can € Nn | ak L66°6 | ce | 
| 0lx621°S [5 _olxezg=2 | beL°O | ©_01x690°6 | sO esos | S2OtxE6E°G. |e Oe eee ee oa 
| Yoveatusn | Soixeie: L | 98L°0 | ‘OLXLON'S | 2 OLxX00h* sa F-OLX0L2°9 o | = 
| 00S | OO€ | O0L | 00S | OO€ | 00L | yySueT una | 
ao = = nn en ne ne ne | nnn nn nnn nn nnn nnn nnn nnn nnn nnn [ncn rnc cnnnccnnn-n| 
. uSTY | MOT | TeuStssestou | 
| 


(SVIGN) uwoyeUTYSY pesetquy 
kq psqyewtysy SAsJZowWedIed syy JO SUOTZETADG Poepueys OBCUdAY 39yuy jo uostueduog Ol th ATGVL 


171 
cases. The standard deviation values for the low and 
high noise/signal ratios were 0.1 and 5.0 for process I 
and 0.01 and 0.1 for processes II and III. 

Computational expense is directly proportional to 
Bhe amount of CPU time used and the two cases 
approximate the low and high ends of the cost of the 
Pours mcaentificatian methods. Table 4.11 presents the 


results of the computational time estimation. 
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4.4 SUMMARIZED RESULTS: 





eeewith (at estimation routines and both 
noise/signal ratio types, increasing the length of the 
input/output sequences always increased the accuracy of 
the parameters estimated when processes I and II were 
tested. 

2) Por==processes = 1 Wand@*rY "with all*e®estimation 
routines, increasing the noise/signal ratio always 
decreased the accuracy of the parameters estimated for 
the same number of data pairs. 

3) Identification using generalized least squares 
did not provide converged process parameter values for 
all «three processes with both noise/signal ratio types 
ano all*’run®* lengths. Convergence was determined by 
CRITE, a user-set comparison variable in UNALIP which 
contained a value of 0.0001 for all tests. 

4) For process III with all estimation routines, 
both noise/signal ratios and all run lengths, estimates 
of the denominator parameters are more reliable than 
estimates of the numerator parameters. 

5) For process III with the generalized least 
squares estimation routine, increasing the number of 
data pairs does not necessarily mean an increase in the 
accuracy of the parameters estimated. This ‘affected arti 
the error measures. Similar behaviour was also noted by 


Isermann [2] and a suggestion that the autoregressive 
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process used to describe the residual sequence is not 
accurate enough was put forth. This accuracy problem 
stems from the truncation of an infinite autoregressive 
process. 

6) For process III with the instrumental variables 
estimation routine and a low noise/signal ratio, the 
error measures were lowest at a run length of 300 data 
pairs. Similar behavior was also noted by Isermann [2] 
except that the standard deviation of the impulse 
response error does’ not ‘follow a steady ‘reduction «in 
error magnitude. Isermann [2] suggested that this 
behaviour results from the model for process III. 

te) For process III with the generalized least 
squares estimation routine and a low noise/signal ratio, 
tne error measures were better than the other three 
estimation routines for a run length of 100 data pairs. 
As the number of data pairs or the noise/signal ratio 
increases, the comparison becomes less and less 
favourable. This could be due to the model for process 
PDI. 

8) With ale estimation routines, both 
Hoase/signal ratio types and all run lengths, the 
standard deviation of the mean squared impulse response 
error related to the mean squared true impulse function 
values for process III exhibits values similar to the 
other two processes even though the parameter estimates 


have large errors. Similar behavior was also noted by 
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Isermann [2]. A possible explanation is that process 
III has six unknown parameters while processes I and II 
have only four. The larger number of parameters would 
allow greater errors in the absolute values of the 
parameter while siida estimating the system 
satisfactorily. 

9) For process III with the generalized least 
squares and the unbiased estimator methods and run 
lengths of 300 and 500, all parameters estimated are 
almost identical, differing only in the. fifth .decimal 
place. No explanation for this-can be given. 

10) The instrumental variable and the unbiased 


estimator methods take the least computational time. 
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4.5 DISCUSSION 


The generalized least squares method showed either 
slow or no convergence of the parameter estimates, thus 
indicating biased parameter estimates. This behaviour 
can be explained by remembering that the infinite 
autoregressive process to describe the residuals must be 
truncated to make it finite for computer computations. 
Thus the residuals are not described accurately enough 
producing a small bias. This appears not to have an 
effect for short run lengths. 

When the noise/signal ratio increases, the 
instrumental variable method and the unbiased estimator 
method produce parameter estimates that converge very 
slowly. Lit smeecounde ber the "result ory paramerver 
approximation step size problems in the two routines. 
The parameter approximation step size determines how far 
the estimation of the parameter should move from the 
previous estimate of the parameter. Many of the 
difficulties encountered when the noise/signal ratio was 
high were due to the step sizes being too small even 
though the estimated parameters were not near the true 
parameters. This could be overcome by increasing the 
complexity of the parameter approximation step size 
reduction program code, but the increased complexity 


must be traded off with an increase in the computational 


time. 
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The maximum likelihood method appears to provide 
the most reliable parameter estimates for different 
numbers of data pairs, noise/signal ratio and process 
type but computational time must also be considered in 
rating the effectiveness of a particualr estimation 
routine. Consideration of the computational time, the 
simulation results show that the preferred 
identification routines would be the unbiased estimator 
and instrumental variable methods, especially when the 


noise/signal ratio is low. 
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CHAPTER FIVE 


PARAMETER ESTIMATION OF MODELS FOR PILOT PLANT UNITS 


5.1 INTRODUCTION 

This chapter contains the work done to estimate the 
coefficients of the pulse transfer functions of two 
Perour plant unitts*% One**unit-. was “a *liquid—-Ti19quid, 
concentric tube heat exchanger. The other was a binary 
distillation column with eight bubble cap trays. 

Aeeoourhn sboth wunits* "exhibit “nonlinear dynamic 
behaviour, the approximate behaviour of both pilot plant 
units can be represented by linear, time invariant, 
difference equations. Therefore the UNiversity of 
ALberta Identification Package (UNALIP) can be applied 
to identify pulse transfer functions that represent 


system behaviour. 
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5.2 DESCRIPTION OF THE PILOT PLANT UNITS 


DISTILLATION COLUMN 

The binary distillation column is 22.86 ecm. in 
diameter with eight bubble cap trays and has glass walls 
for viewing. Dra yrs ta Cringe ets =-3 0-10" oma All 
temperatures, flow rates, and compositions measured on 
the process are transmitted to an IBM 1800 digital 
computer and logged for data acquisition and/or control 
action by means of the standard Direct Digital Control 
(DDC) program [26]. The distillation column is equipped 
with analog controllers to control flow rates, reboiler 
and condenser eee and the Bee and reflux 
temperatures. The Plowscontroller e6t  pointsiare 
manipulated by means of the DDC progran. A detailed 
description Oofmathe distil Wation column and the 
associated equipmenet is given by Svrceek [22] and Pacey 
i723). Figure 5-1 is” a simple schematic ‘of the 


distillation column configuration. 
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FIGURE 5-1 
SIMPLE SCHEMATIC OF THE DISTILLATION COLUMN 
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HEAT EXCHANGER 

The double pipe heat exchanger is operated in 
countercurrent manner with both streams being water. 
All temperatures and flow rates measured on the process 
are accumulated by the Direct Digital Control programs 
on the IBM 1800 digital computer, except that the tube 
Side flow rate and the outlet temperatures for the tube 
Side and shell side fluids are accumulated by the High 
Speed Data Acquistion (VCDAQ) package [24]. Tube side 
flow is regulated by a pneumatic control valve with the 
pneumatic signal manipulated by the IBM 1800 digital 
computer executing the heat exchanger dynamic testing 
progran. A detailed description of the heat exchanger 
and associated equipment is given by Wildman [24] and 
Lees [25]. Figure 5-2 is a simple schematic of the heat 


exchanger configuration. 
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SCHEMATIC DIAGRAM OF HEAT EXCHANGER 
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5.3 EXPERIMENTAL PROCEDURE 


DISTILLATION COLUMN 





The object of the experiment on the distillation 
column was to obtain the four pulse transfer functions 
relating the inputs, the steam and reflux flow ratesshto 
the outputs, the top and bottom compositions. The 
transfer functions would then be used in a control study 
being conducted by Bilec [26]. The steam and reflux 
flow rates were varied, using a pseudo random binary 
sequence as the forcing signal to the flow control 
valve. Pseudo random binary sequences of length 63, 
double 31 and 127 were used. The flow rates and the top 
and bottom compositions were logged via the DDC program 
er tthe ¢xIBM ! 9 1:800:- The resulting input/output data 
sequences were analyzed and the pulse transfer function 
coefficients determined by the identification routines 
RerkMei3.10}, IVOFF {3.11} and UBIAS {3.12}. A detailed 
description of the experimental procedure is given by 


Bilec [26]. 


HEAT EXCHANGER 

The object of the experiment on the heat exchanger 
was to obtain the pulse transfer function relating the 
input, the tube side flow rate, to .the output, unersnell 
side fluid outlet temperature. The tube side flow rate 


Was varied, using a pseudo random binary sequence as the 
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forcing signal to the pneumatic control valve. The tube 
side flow rate and shell side fluid outlet temperature 
were logged by the VCDAQ program for subsequent 
analysis. The data was collected over 50 seconds at 10 
points per second. Pseudo random binary sequences of 
pengense 255 piwithsca clock intervali;ofurl ,aandesizuwithea 
elock interval of 10, were used. The clock interval is 
the number of sampling periods that the pseudo random 
binary sequence forcing signal to the pneumatic control 
Valve »~Maintains -a& certain state’. [nheseresulvlansg 
input/output data sequences were analyzed to determine 
une pudse transfer functions using | each} of \othe 
fdentitication techniques. STIME {3.9}i, “ASTRM? {3.10}, 
PVOFSei seit and “UBIAS {3.12}. A detailed description 


of the experimental procedure is given in Appendix A. 
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5.4 RESULTS 


DISTILLATION COLUMN 





Tables 5.1,to 5.4 are tabulations of the system 
parameters’ efor «thei “identification runs involved. 
Appendix C has the plots of the input and output data 
fer seach test. The averaged pulse transfer functions 


obtained for the distillation column were 
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where: 1) Xp is the top composition (weight percent 
methanol) 
2) X, is the bottom composition (weight percent 
methanol) 
Sole ts the steam “flow ratep to 1the~ reboller 
(g/s) 
4) RE is the reflux flow rate (g/s) 

Comparison of the gains of the estimated transfer 
functions with those established by step testing showed 
that there was large disagreement. Consequently in his 
control study, Bilec [26] used the gains from the open 
loop experimental tests and the time constants 
established by this identification study. The transfer 


functions that he used were: 











a SpE aie eepre 
X,(s) 1037.4s + 1 1AM. See Slush 
; Chiee 
Xy(s) -~0.6310 Pee see RE(s) 
859.3s + 1 Pio. | sok a 


After SapudyLug «Vine problem of poor gain 
realization, it was noticed that the transfer functions 
with the poorest gains were the transfer functions that 
had the largest time constants. A simulated system was 
then studied to analyze the results from the different 


identification techniques to ascertain the cause of poor 
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estimation of the transfer function gains. The 


Simulated system was 


Se ioe ene OSs (5.4) 


fe. 100s 1 = 029 9ra 





The actual gain of the simulated system is 1. The input 
of the simulated system was a pseudo random binary 
sequence of length 31. The clock interval used was 1. 
miemvOtal number of <input ‘points was 500. The output 
sequence was corrupted with white noise that had zero 
mean and a standard deviation of 0.1. The results of 


the experiment are given in Table 5.5. 


TABLE 5.5 Identified Gains of Simulated System 
Using Each Estimation Technique 


| estimation | gain | 
| routine | | 
| ------------|---=------1-| 
|  STIME [Pe0e85 x.10.ueun| 
ASTRM 0.118x10, | 
IVOFF 0.155x10°, | 


| UBIAS | 0.154x10- 
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The results in Table 5.5 confirmed the suspicion 
that there were no errors in the implementation of the 
identification techniques but that the difficulty was 
likely the particular choice of the pseudo random binary 
sequence input used for the tests. Consequently, it was 
decided to assess the effect of a change in the clock 
interval of the input pseudo random binary sequence upon 
the estimated gains. The clock interval was varied from 
5 to 20 in predetermined stages to see if this would 
improve the gain estimates [2]. The length of the 
pseudo random binary sequence was maintained constant at 
a length of 31. The total number of data pairs was 
again selected at 500 and the calculated system output 
corrupted with zero mean and 0.1 standard deviation 
white noise. The results from identification performed 
arose resulting data, pairs for ,particular “clock 
intervals are given in Table 5.6. 
Table 5.6 indicates, that increasing the clock 
interval from 5 to 15 increases the accuracy of the 
estimated gains but further increases appear to have a 


negative influence. 
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TABLE 5.6 Effect of Clock Interval on the Estimated 
Gain of a Simulated System 


| clock | estimation | gain | 
| interval | routine | | 
| ----------| == ---------- fe lortec i tant aie = 
| | STIME \ie0 146 xX1 Oneene 
| 5 | ASTRM | 1.038 | 
| | IVOFF | 0.501 | 
| | UBIAS [Ook | 
| ---------- | ------------ | ---------- r-| 
| | STIME | 0.494x10 | 
| 10 | ASTRM \oreO 36 | 
| | IVOFF Li eteaOal 0 | 
| | UBIAS [02957 | 
lana ee ee ee ie ae ol 
| | STIME | 1.100 | 
| 15 | ASTRM bates | 
| | IVOFF he Dee | 
| | UBIAS | 1.003 | 
SSsseesesons | ------------ ~Saaeorecaoss 
| | STIME iy41.355 | 
| 20 | ASTRM | 0.954 | 
| | IVOFF (po. 862 | 
| | UBIAS | 0.862 | 
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HEAT EXCHANGER 


The pulse transfer functions obtained when a pseudo 
random binary sequence of length 255 and clock interval 


of 1 sampling interval (0.1 sec.) was used for the input 


forcing function were: 
1) Generalized Least Squares (STIME {3.9}) 


To(z) -0.00128927 || 


= (5.5) 
cz) 0.992727 


where: 1) Fo is the tube side flow rate (liter/min) 
a8 Ty is the shell side fluid outlet temperature 


(degrees C) 


2) Maximum Likelihood (ASTRM ec 


Tet ey =puoog0522% | rey 


Ree 1 0.996327! 
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3) Instrumental Variables (IVOFF {3.11}) 





a 8A) 6 oqousest? ere 
Po(2z) 4. 9.979727! 

4) Unbiased Estimator (UBIAS {3.12}) 
T,¢2) Peo .ciosgzm = Ey 
Foz) 4 _ 9. 980227! 


The input forcing function values are generated between 
Zeroe ang plus or minus one. These values are then 
converted to controller output values acceptable to the 
control value. The initial steady state tube flow rate 
ese este l/min. The input forcing function would cause 
the tube flow rate to vary + 11.35 l/min. depending on 
the sign of the input forcing function. 

The pulse transfer functions obtained when a pseudo 
random binary sequence of length 31 and clock interval 
of 10 sampling intervals (1.0 sec.) was used for the 


input forcing function were: 
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1) Generalized Least Squares (STIME {3.9}) 


B. GZ) = 

Ss -0.009326 
oe SE (5.9) 
Fp(z) | SARE! 


2) Maximum Likelihood (ASTRM ttt Obs) 


ea 2) -16 

S -0.018 
aie 55z (5.10) 
F(z) ors 0.984327! 


3) Instrumental Variables (IVOFF {3.11}) 


(Need GEAR -16 
S i -0.01855z (5.11) 


Fo(Z) 4 _ 9.984327! 
4) Unbiased Estimators (UBIAS) 


eh al OA) -16 
ee (5.12) 


ay 0. 984 bz5 | 


The input forcing function values are generated between 
zero and plus or minus one. These values are then 
converted to controller output values acceptable to the 


control valve. The initial steady state tube flow rate 
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Was ec.) 27min. The input forcing function would cause 
the tube flow rate to vary + 11.35 l/min. depending on 
the sign of the input forcing function. 

Appendix B contains the plots of the experimental 
input and output data for each test run. The steady 
State energy balance, initial test conditions and a 
Summary of the psuedo random binary sequence conditions 
are given in Appendix A. 

Wildman [24] determined continuous’ transfer 
fmunctwvons by "pubse "(esting using *four different Gnput 
pulse shapes. Averaging of the time constants (time 


delays) and gain values gives the transfer function 


1.14s 


Ty (s) 2 Sieger 


(5.13) 
Fs) hes oe Pet a 


Maresthe Sampling time of 0.05 seconds, the pulse 


transfer function is 


Ty (2) x 0.016742 || acing 





Bee Z)e 7 _ 9.984325. 


As can be seen, the agreement between the average 
transfer function representation established by Wildman 
and the transfer functions estimated by the ASTRM 


meeioy, IVOFF {3.11} and UBIAS {3.12} routines are good 
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for both experimental runs. However the best estimates 
resulted using the psuedo random binary sequence of 
length 31 and clock interval of 10 sampling intervals 
(1.0 sec.). Excellent gain extimates were obtained. 

It is to be noted that the generalized least 
squares estimated model was not consistent with the 
models established using the other three identification 
routines nor with the model Wildman determined [24]. 
Similar behavior was noted in Chapter Four. 

These results are thus consistent with the model 
estimates in Chapter 4 established using the different 


routines. 
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5.5 DISCUSSION 


The estimation of the pulse transfer functions by 
the UNiversity of ALberta Identification Baan ed 
(UNALIP), in general, was acceptable. The generalized 
least squares indentification (STIME {3.9}) routine was 
consistently poorer than the other three routines, 
producing a time delay and gain values that were 
unacceptable. This is the reason that only the maximum 
likelihood (ASTRM {3.10}), instrumental variables (IVOFF 
{3.11}) and unbiased estimator (UBIAS {3.12}) routines 
were used to estimate pulse transfer functions for the 
distillation column. 

For specific transfer functions, the estimation 
routines for individual disturbance types produced 
results that wereilivery ‘close’s'to each other: For 
different disturbance types, the results were not always 
consistent. Pseudo random binary sequence input signal 
lengths of 63 and double 31 produced comparable results 
but the estimates from the pseudo random binary sequence 
input signal of length of 127 were not comparable to the 
other tests. On the basis of the simulation results in 
Chapter Four it was concluded that the results from the 
pseudo random binary sequence length of 127 should be 
used because generally parameter estimate accuracy 


increased if the number of data pairs increased. 
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A possible reason for the poor gain estimation for 
the distillation column is that for systems with large 
time constants, the clock interval must be large so that 
the process is able to ‘'see' the input disturbance. A 
criterion for a good input sequence is that it excites 
all possible states of the process [27]. This was 
probably not the case for the distillation column tests. 
Unfortunately runs to test this were not done on the 
distillation column but a type of check was accomplished 
on the heat exchanger. A pseudo random binary sequence 
of length 31 and clock interval of 10 gave better 
results that the pseudo random binary sequence of length 
255 and clock interval of 1. Also, the estimates of the 
heat exchanger transfer function were very similar for 
the estimation routines ASTRM {3.10}, IVOFF {3.11} and 
UBIAS {3.12} which is consistent with the previous 
Simulation results. 

There was good agreement between the estimated time 
delays for ®the “distillation column pulse transfer 
functions and the time delays determined by Bilec [26] 
from step testing. This Wwastinotisthe case Torsther near 
exchanger but from the range in the values of the time 
constants determined by Wildman [24], it is reasonable 
to conclude that the estimation of small time delays is 
difficult for a real process. Wildman [24] estimated 


the time delay to be equal to 11 sample intervals or 
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0.55 seconds, which is close to the 16 sample intervals 
or 1.60 second time delay estimated in this study. 
It can be concluded that when estimating the pulse 
transfer function of an actual process that: 
1) as many tests as possible with different input 
sequences should be performed using 
a) different clock intervals 
b) different length of input sequences 
c) different types of input sequences 
d) different magnitude of input sequence 
changes from steady state 
2) as many identification methods as possible 
Should be applied to the data 
3) in deciding upon the final estimated pulse 
transfer function from the functions provided 
by the identification routines, care should be 
taken to compare the estimates from various 
routines and use only those estimates that 
experience dictates «are; reasonable for the 
process under consideration to obtain an 
average transfer function. 
MoecpendaxesD 1s ‘a collection) of jguideliinesastor 
performing identification of an unknown system. These 


guidelines were established during the work on the pilot 


Paanthunits:: 
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CHAPTER SIX 


CONCLUSIONS AND RECOMMENDATIONS 


The following conclusions can be drawn from this 


study: 


1) 


2) 


3) 


4) 


On the basis of analysis of the experimental tests 
and generated data sequences, the preferred 
estimation routines were the instrumental variable 
and the unbiased estimator methods, considering both 


parameter estimation accuracy and computer cost. 


The maximum likelihood technique, provided the most 
accurate parameter estimation for the tests using 


generated input/output system values. 


The instrumental variable and the unbiased estimator 
methods provided the best parameter estimates using 


the experimental data from both pilot planLeuni ves 


Increasing the number of input/output data pairs used 
by the estimation routines does not always increase 
the accuracy of the parameter estimates. The order 

of the system being studied has an effect on the 


accuracy of the parmeter estimates. 
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Increasing the noise/signal ratio for the three 
Simulated processes decreased the accuracy of the 


parameter estimate for all estimation routines. 


6) The accuracy of the parameter estimates decreases for 


i) 


all estimation routines as the order of the generated 


System increases. 


In a first order system, increasing the clock 
interval of the input sequence increases the accuracy 
of the estimated gain if the time constant is large 
compared to the sampling interval. The -clock 
interval must be varied until gain estimates are 
comparable to open loop values of the gain. Tt ss 
possible to increase the clock interval of the input 
Signal too much thus giving poor estimates of the 


gains. 


The following are some recommendations for future 


work: 


1) The UNiversity of ALberta Identification Package is 


now implemented on the University of Albert's Amdahl 
470V/6 digital computer. The package should be 
implemented on the Department of Chemical 
Engineering's HP1000 distributed network digital 


computer. This would allow easy access Lo tne 
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2) 


3) 


4) 


<p) 
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package for engineering studies conducted in the 
department. One disadvantage would be the reduced 


response time obtained by the package on the HP1000. 


All the estimation routines in UNALIP are of the 
off-line type, an improvement in the arrangement 
would be to convert the routines to an on-line type 
for implementation on the HP1000 to allow real time 


model estimations to be included in control studies. 


All the estimation routines can be converted to 
recursive forms which would save computational time 
and save core memory space. This would make 
conversion of UNALIP from the Amdahl mainframe to the 


HP minicomputer more feasible. 


Some theoretical work in generalized least squares 
could be done on the problem of convergence of the 
parameter estimates to the neighbourhood of their 
true values [8]. The generalized least squares 
routine could then either be changed to promote 
convergence or eliminated as an estimation routine vin 
UNALIP since it is currently doing an unsatisfactory 


HOD. 


More testing of the identification package should be 


carried out, especially the use of high pass filters 
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and differencing used in the routine DMOVE. The 
testing would determine if these options help in the 
estimation process and if ° so, how do they aid 


improved estimation. 


When the identification package is implemented on the 
HP1000 digital distributed digital computer network, 
it is recommended that a time series analysis program 
be included since there is no program of this type 


currently on the system. 


An@vinvestigatvion» #into.sthe -effect’Yof the clock 
interval on the estimation process should be carried 
out. This study could look into: 
1) what value of the clock interval will give the 
best estimates for different processes 
2) does the clock interval affect the estimation of 
the process's time delay 


3) when should the clock interval be altered 
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This appendix contains the detailed description of 

the experimental procedure to obtain the data from the 
heat exchanger. This data was used by the estimation 
routines in UNALIP to provide the results discussed in 


Chapter Five. 
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Experimental Data for Test A 


Forcing Function - pseudo random binary sequence 


- 255 cycles/pulse width 
- clock interval equal to one 
Sampling interval 


Response Function - shell liquid outlet temperature 


Flow Conditions - countercurrent, water-water 
Heat Balance Data 


mooetluid flow (]1/min) - 
Shell fluid flow (1/min) - 


Tube Fluid inlet temperature (degrees C) - 
Tube fluid outlet temperature (degrees C) - 


Shell fluid inlet temperature (degrees C) - 
Shell fluid outlet temperature (degrees C) - 


Experimental Data 


Duration of forcing function (seconds) - 
Total time for data accumulation (seconds) - 
Number of data points accumulated per second 


Amplitude of the forcing function (percent) - 


cor 
45. 


Wc 
20r. 


tds 
68. 


50 
50 


10 
5 


214 


on) 


“a Wo W 





|) . a, a, 
a ~ o Phat 7 


‘ a ¥ A 
be v 
2 
' 
: 










pome ny rhe ~9an iit er ry 
ata < 

cg Oh aap a bavraaas A905; 
i, yV'S ose} * watt _— " i 


o>" ie 
(rad taliuo Hide ht if "got 






























as ‘oe 
Dies oe ire 

w«IoJeu , Teenthe*e4 avog = _ 
Vy i eee 

(pdents wofl et cy oe 

CRLELET BOSE b tid ae Ne 


3 £6o6"g90) stags asenet ‘deter baud % 
ra (i) Boetes e) yt se4oqnaa ah a] a 


J \ee dee rere eqs. total bie tea 
. 3 ) sess egeed detsvo lace: — 


“— (ebproctee) adidony? gatos 16.96) 
. (ebeosse) aol’ eivnysos goed tot ote 
ee +c : Sa JS LUBAWSUS etaiog sind. 7 


roisoou? paLore ont to 





a he 
Experimental Data for Test B 
Forcing Function - pseudo random binary sequence 
- 31 cycles/pulse width 
- clock interval equal to ten 
Sampling intervals 
Response Function - shell liquid outlet temperature 


Flow Conditions - countercurrent, water-water 


Heat Balance Data 


mube: fluid flow (1/min) mee iT 
Shell fluid flow (l/min) - 45.4 
Tube Fluid inlet temperature (degrees C) - 12.3 
Tube fluid outlet temperature (degrees C) mio ea eh 
Shell fluid inlet temperature (degrees C) Oi | 
Shell fluid outlet temperature (degrees C) - 69.0 
Experimental Data 

Duration of forcing function (seconds) - 50 
Total time for data accumulation (seconds) - 50 
Number of data points accumulated per second 

- 10 


Amplitude of the forcing function (percent) - 5 
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Appendix —B contains = the plots of =“the=-tnput and 
puepuce gata for -each-test_.conducted Jon. the heat 


exchanger. 
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Tube Flow Rate Response for Test A 


Figure B-1 
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Tube Flow Rate Response for Test B 


Figure B-2 
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Shell Outlet Temperature Response for Test A 


Figure B-3 
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Shell Outlet Temperature Response for Test B 


Figure B-4 
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Appendix C contains the plots of the input and 
output data for each test conducted on the distillation 


eolumn. 
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Reflux Flow Rate Response for PRBS 


Figure C-1 
Disturbance of Length 63° 
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Reflux Flow Rate Response for PRBS 


Figure C-2 
Disturbance of Length Double 31 
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Steam Flow Rate Response for PRBS 


Figure C-3 
Disturbance of Length 63 
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Figure C-4 Steam Flow Rate Response for PRBS 
Disturbance of Length Double 31 
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Bottom Composition Response for PRBS 


Figure C-5 
Disturbance of Length 63 using Steam 
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Bottom Composition Response for PRBS 


Figure C-6 , 
s Disturbance of Length 63 using Reflux 
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Bottom Composition Response for PRBS 


Figure C-7 , 
Disturbance of Length Double 31 using Steam 
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Bottom Composition Response for PRBS 


Figure C-8 : 
Disturbance of Length Double 31 using Reflux 
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Top Composition Response for PRBS 


Figure C-9 r 
Disturbance of Length 63 using Steam 
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Figure C-10 Top Composition Response for PRBS 
Disturbance of Length 63 using Reflux 
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Figure C-11 Top Composition Response for PRBS 
Disturbance of Length Double 31 using 


Steam 
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Figure C-12 Top Composition Response for PRBS 
Disturbance of Length Double 31 using 
Reflux 
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Figure C-13 Reflux Flow Rate Response for PRBS 
Disturbance of Length 127 
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Figure C-14 Steam Flow Rate Response for PRBS 
Disturbance of Length 127 
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Figure C-15 Bottom Composition Response for PRBS 
Disturbance of Length 127 using Steam 
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Bottom Composition Response for PRBS 


Figure C-16 
Disturbance of Length 127 using Reflux 
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Figure C-17 Top Composition Response for PRBS 
Disturbance of Length 127 using Steam 
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Top Composition Reponse for PRBS 


Figure C-18 
Disturbance of Length 127 using Reflux 
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GUIDELINES FOR PERFORMING PSEUDO RANDON BINARY SEQUENCE 


IDENTIFICATION 


This appendix contains a collection of guidelines 
deveronoed. during, the. testing of _ thes two “pilot. wuniits 
investigated during this thesis project that made 
identification of process units more systematic and less 
time consuming. The guidelines pertain only to open 
Troops vSTems Since .tonly, this whype,.ol syst emwvas 


investigated. 
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A) INPUT SIGNAL CHARACTERSITICS: 

1) amplitude of input signal - this factor can be 
established by running open loop tests to determine 
system output response to step changes, up and down, 
from normal input conditions. Ihew ODJEGLLVeY 1s 4vO0 
produce an output response that covers as much of the 
high and low ends of the spectrum of responses as 
possible without running into physical limitations on 
the system or operating constraints. UNALIP assumes the 
systems to be ‘identified are linear. Non-linear 
behaviour becomes more apparent as the extremes of 
system response are reached and hence the validity of 
the linearity assumption can be checked more readily. 

2) length of input signal - the length of the input 
signal depends on the dynamic behaviour of the system 
under consideration and the sampling interval. An 
observation made during the testing of the distillation 
column was that the length of the input signal should be 
at least as long as the largest time constant of the 
system under eBiatdcratton. The’ step response of a 
first order system is only 63% complete after one time 
constant. Thus shortening the length of the input 
signal would omit a significant amount of response 
information. Sometimes system considerations or 
computer limitations do not make this possible. If this 
is the case then the general rule is that the longer the 


input signal length, the better the estimation. Some 
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work on sampling theory is discussed by Wozencraft and 
Jacobs [10] and can be applied in this area but was not 
done in this thesis. 

3) sampling interval - the sampling interval is 
crucial to obtaining good system estimates and depends 
on system dynamics and physical limitations of the 
sampling mechanism. The sampling interval should be 
one-tenth as long as the smallest time constant of the 
system. When in doubt, the analysis can be repeated 
WMoecneescverai~ Urial= sampling=eintervals**as"a- simple 
sensitivity test for the length of sampling interval. 
In’ the® choice” of “sampling interval “it is” the* output 
foresce=onave’Ys important’,-* and §at's “varvVance @*s roird 
approach a minimum value as the interval is shortened. 
Another guide for sample interval values is experiments 
done on similar equipment, as was the case with the heat 
exchanger experiments carried out in this thesis. The 
tradeoff of input sequence length vs sampling interval 
COccours sat! the "total ‘number’ of? points “is "fixed ‘ory if 
large differences in system time constants exist. Thus 
the lower limit of a sampling interval for a system 
depends on system dynamics or computer storage size. 

4) clock interval - the clock interval of a psuedo 
random binary sequence will depend on the size "of "the 
sampling interval and the magnitude of the system time 
constants. If the sampling interval is small and some 


of the time constants are large, then the clock interval 
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will have to be large to allow the system to 'see' the 
input disturbances. This will improve the estimates of 


the system gain which can be checked by doing open loop 


gain estimations. 


B) SYSTEM CHARACTERISTICS: 

1) model order - an initial estimate of the model 
order of a system can be determined from the 
characteristic shapes of the open loop response of the 
System’ to» step changes’ "in the” input Signal: Tess 
carried out on similar systems can also be used as a 
Starting point. UNALIP offers a test to determine which 
model order>*is correct. If there is no dependence 
between the prediction errors and the input, then the 
assumed model order is adequate. The best model is 
determined by comparing the mean square error of 
different model orders. If the higher order model is 
significantly better than the lower order model, the F- 
ratio test will indicate this. If the covariance matrix 
has a determinant of zero then the model order is too 
high. For good estimates of the system model order, it 
is best to use several estimation routines to check on 
the results as was observed in the tests on the heat 
exchanger. 

2) time delays - as with model order, initial 


estimates of the time delay can be obtained from open 
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loop experiments and experience with similar systems. 
UNAL1IP also has tests to determine if the time delay is 
correct. If there is no dependence between the residuals 
and the input then the time delay is correct. For time 
delays it is best to use more than one estimation 
routine to check the results, as was observed in the 


tests on the heat exchanger. 


C) PROBLEMS: 

1) noise - the higher the noise/signal ratio, the 
more difficult estimation becomes. Tnew order sO) waa 
system determines the amount of noise that can be part 
of a system's output before parameter estimation 
problems occur. The simulated first order system, with 
a noise/signal ratio of 5.0 had similar estimation 
success to a simulated second order system with a 
noise/signal ratio of 0.1. The higher the noise/signal 
ratio, the higher the failure of the estimation routines 
to converge or agree with other estimation routines. 

2) model order - the larger the model order, the 
more difficult estimation becomes. Anything larger than 
a second order system significantly increased the 
calculation time required and resulted in models that 
gave good estimates of steady state gain but poor 


estimates for individual parameters. 
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3) time delays - the smaller the sampling interval, 

the more difficult time delay estimation becomes. This 
problem was noted in the tests done on the heat 
exchanger. This point should be remembered when placing 
confidence in time delay estimation of systems with fast 
dynamics, although as time delays become smaller their 


relative importance probably decreases. 
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